Reminder: Day 2 Schedule

9:05-10:20 - Session 1

  • Day 2 Tutorial, Part 1

10:20-10:30 - Snack Time

10:30-11:55 - Session 2

  • Day 2 Tutorial, Part 1 continued

12:00-12:45pm - Lunch @ Pavilion

1:00-3:15 - Session 3

  • 1:00 - 3:00pm Day 2 Tutorial, Parts 2-3
  • 3:00 - 3:15pm Day 2 Tutorial. Part 4, save your work

Part 1. Mathematical modeling

Run the package installation script again and load the packages you need for this morning.

source("installscript.R")
require(geodata);require(dismo);require(maps);require(here)
## Loading required package: geodata
## Loading required package: terra
## terra 1.8.54
## Loading required package: dismo
## Loading required package: raster
## Loading required package: sp
## Loading required package: maps
## Loading required package: here
## here() starts at /Users/jblois/Documents/GitHub/biodata_shortcourse/development

Now, make sure your workspace from yesterday is loaded. There are several options for doing this: - go to the Documents folder, choose the .Rdata file that you saved yesterday afternoon, and open it in RStudio. It should look exactly like it did yesterday when you left it. - Open RStudio, go to ‘Session’, then “Load Workspace”. Navigate to your saved workspace and open it. - Load it directly within RStudio using the example code below.

load(file = "Blois_Day1.RData") # Replace the "Blois_Day1" with your filename. 

Introduction to species distributions

Knowledge about where a species occurs, and where it doesn’t, is fundamental to conservation biology. Over the decades, scientists have collected a lot of different species occurrences - observations of a species at a particular place and time. But, because we can’t be everywhere all the time, there are many dots on the map where a species is likely to occur, we just haven’t observed it there yet. One way to predict where a species can be found is to build a model of the species distribution (also called the species range).

Today, we’re going to use the climate data you downloaded and the species occurrence data you constructed yesterday to model (determine) the influence of the climate variables on your species’ range.

Predicting species presence from climate variables

Species distribution models describe the climate where a species currently lives (its climate “envelope”) based on the values at all its known occurrences, and then use that to predict the suitability of all the places it has not been observed at. We’ll try two different ways of making a species distribution model.

First, take a look at your environment to see if you have the Data object climate loaded from yesterday. If so, skip the next chunk of code. If you need to rebuild the climate raster stack, run the following lines:

climate <- worldclim_global(var = 'bio',res = 2.5,path = here())
climate <- stack(climate)
names(climate) <- unlist(sapply(1:19,function(x) paste0("bio",x)))

Run the following line to make the climate object smaller:

climate <- rast(climate)

Bioclim model

We can use the function bioclim to make a simple species distribution model known as a “climate envelope model”. The Bioclim model, specifically, compares the values of the climate variables at each point on the map to the values at the locations where the species is known to occur.

Reminders: - the object df stores the occurrence points you collected yesterday. - To see what each of the 19 bioclimatic variables means, look at https://www.worldclim.org/data/bioclim.html. Temperature measurements are given in tenths of a degree Celsius; precipitation is in millimeters.

m1 <- bioclim(stack(climate),df[,c("lon","lat")])

Now that we have the relationship between climate variables and your species known occurrences, we can predict its likelihood of occurring for every point on the map using the predict function.

p <- predict(climate,m1,progress = "text") # this line of code takes a minute or two to run.
## |---------|---------|---------|---------|=========================================                                          

We will save and reload this file to get it out of RAM:

p <- reload_raster(p)

We can plot the predicted ‘climate suitability’ (a prediction of the species distribution) just like we plotted the bioclim environmental layers:

plot(p, main = "Bioclim model", legend.args = list(text = 'Suitability', side = 4, 
         font = 2, line = 2.5, cex = 0.8)) 

The suitability score ranges from 0 to 1. A point that had a suitability of 1 would have all its climate values exactly at the median values for the occurrence points, so that doesn’t usually happen; the highest values in this model are around 0.6.

This model seems to think that the most suitable habitat for blue morpho butterflies is all the areas in the world with tropical rainforest. Let’s zoom the map in on the range containing the actual data points, like we did before…

xrange <- range(df$lon,na.rm = T) + c(-1,1)
yrange <- range(df$lat,na.rm = T) + c(-5,5)

plot(p, main = "Bioclim model prediction",
  xlim = xrange, ylim = yrange,
  legend.args = list(text = 'Suitability', side = 4, 
         font = 2, line = 2.5, cex = 0.8)) 

…and also put in the occurrence points.

plot(p, main = "Bioclim model prediction with occurrences",
  xlim = range(df$lon,na.rm = T), 
  ylim = range(df$lat,na.rm = T),
  legend.args = list(text = 'Suitability', side = 4, 
         font = 2, line = 2.5, cex = 0.8)) 
points(df$lon,df$lat,col = rgb(1,0,0,0.5),pch = 20,cex = 0.5) #make the points smaller so they don't get in the way too much

The Bioclim prediction doesn’t appear to do a particularly good job of predicting the species distribution for the blue morphos. Notice how some of the most suitable areas don’t have any specimens at all, and many of the specimens are from places with low modeled suitability. Either climate isn’t the main factor in where the blue morphos can live, or the sampling is not very good, or both.


Question: How well did the Bioclim model work for your species? Save the plots to your folder.


Linear model

There are other ways to make a species distribution model than the function bioclim. One way is to develop a “linear model”. The equation for a line (y = mx + b) is a linear model. In terms of species distributions, a linear model is an equation that describes how the value of each of the climate variables in a particular spot on the map influences the likelihood that the species will occur in that spot.

To develop a linear model we’re going to need to make some points in places where your species isn’t present (or at least has not been documented to occur), so that we can compare them to the places where the species was actually found. The other model just used the whole background supplied by the maps of climate.

First, get the values of the climate variables at the location of each of the known species occurrences with the function extract from the raster package. This function takes two arguments, a RasterStack (ie, the object climate) and a set of coordinates (GPS points, ie, the object df), and returns a matrix of all the climate variables at all the locations.

presvals <- extract(climate, df[,c("lon","lat")])
presvals <- presvals[,-which(names(presvals) == "ID")]
dim(presvals)
## [1] 561  19

Now we need to make up some long/lat points for absence data – we’ll try to use approximately use 1000 or the number of presence points, whichever is greater (though if you have more than 10,000 presence points, don’t use more than that). For the blue morphos, there are 561 presence points, so I will develop a dataset with 1000 absence points.

Here, we are “masking” out (not using) any point that has a known occurrence in it, and only choosing other points.

set.seed(0) #initialize random number generator
backgr  <-  randomPoints(mask = raster(climate$bio1), n = max(nrow(df),1000), 
                         p = df[,c("lon","lat")],
                         ext = extent(rbind(range(df$lon),range(df$lat)))) 

You can map the absence points we just made up.

map("world2",col = "darkgray",
     xlim = range(transform_lon(df$lon),na.rm = T) + c(-1,1), 
     ylim = range(df$lat,na.rm = T) + c(-1,1))
map.axes()
points(transform_lon(backgr[,"x"]),backgr[,"y"],col = "lightblue", pch = 20) 
#the absence points cover the whole space
points(transform_lon(df[,"lon"]),df[,"lat"],col = "red",pch = 20)

The blue circles are the “potential absence” long/lat points that we just made up for locations where the species is not documented to occur. We’ll extract the climate values at these “absence” points and assemble them into a table with the “presence” values.

absvals <- extract(climate, backgr) #extract climate values at those points
pb <- c(rep(1, nrow(presvals)), rep(0, nrow(absvals))) #make a vector that indicates whether the species is present at that point or not
sdmdata <- data.frame(cbind(pb, rbind(presvals, absvals))) #bind the presence and absence values together into a single data frame
head(sdmdata) #rows with pb = 1 means the species was found at that location
##   pb     bio1      bio2     bio3      bio4   bio5   bio6   bio7     bio8
## 1  1 21.86033 10.082000 58.83520 230.33929 30.172 13.036 17.136 24.50400
## 2  1 25.95267 10.109333 76.81864  77.31625 32.812 19.652 13.160 25.28133
## 3  1 22.46550  9.664333 89.51773  39.99062 27.904 17.108 10.796 22.17667
## 4  1 25.16517 11.985667 70.30541  48.31260 33.444 16.396 17.048 25.26467
## 5  1 25.08783 10.650333 71.02116 104.88419 31.860 16.864 14.996 25.78600
## 6  1 23.55300  9.641333 62.13801 191.79256 31.060 15.544 15.516 24.98867
##       bio9    bio10    bio11 bio12 bio13 bio14    bio15 bio16 bio17 bio18 bio19
## 1 19.00733 24.72067 19.00733  1292   199    25 63.85461   580    90   566    90
## 2 26.07533 27.06933 25.11600  2853   435    25 54.75378  1106   176   562   941
## 3 22.71400 22.77600 21.86400  4054   464   244 23.67735  1338   749   858  1154
## 4 24.68200 25.59667 24.46933  2130   326     8 69.36523   953    54   301   128
## 5 23.61200 26.00333 23.57533  1975   291    46 58.45269   844   152   640   186
## 6 21.63600 25.83333 21.09533  1192   179    40 50.85691   514   146   358   146
tail(sdmdata) #rows with pb = 0 means there's no occurrence at that location
##      pb     bio1      bio2     bio3      bio4   bio5   bio6   bio7     bio8
## 1556  0 26.47483  9.229667 85.14452  26.48380 32.040 21.200 10.840 26.43600
## 1557  0 25.88550  8.684334 86.11993  49.58579 30.860 20.776 10.084 26.00333
## 1558  0 24.04933 12.374667 69.56750  61.73869 32.324 14.536 17.788 24.16267
## 1559  0 22.28983 13.867666 69.00710 157.46184 30.800 10.704 20.096 23.34600
## 1560  0 20.04800 11.283999 88.29421  26.40729 26.448 13.668 12.780 19.81400
## 1561  0 22.65033 12.238000 69.80379 137.51379 30.204 12.672 17.532 23.75933
##          bio9    bio10    bio11 bio12 bio13 bio14    bio15 bio16 bio17 bio18
## 1556 26.53667 26.80333 26.15000  2309   295    89 38.87713   837   297   478
## 1557 25.35800 26.28600 25.14667  2995   324   210 15.09136   928   669   769
## 1558 23.36067 24.60333 23.14267  1915   335     6 77.64484   930    43   426
## 1559 19.98133 23.43333 19.98133  1740   299     8 76.57098   829    38   684
## 1560 19.91067 20.29933 19.75067  2956   414    99 42.10091  1129   355   733
## 1561 20.66333 23.75933 20.66333  1526   240    14 69.13038   687    56   687
##      bio19
## 1556   405
## 1557   688
## 1558    77
## 1559    38
## 1560  1068
## 1561    56

These presence and absence data can now be used to make a model using the function glm() (which is short for ‘generalized linear model’).

m2 <- glm(pb ~ bio1 + bio2 + bio3 + bio4 + bio5 + bio6 + bio7 + bio8 + bio9 + bio10 + bio11 + bio12 + bio13 + bio14 + bio15 + bio16 + bio17 + bio18 + bio19, data = sdmdata)

This is a linear model, so you can look at the p-values to assess whether each variable has a significant influence on the species distribution. (Remember that you can check what the different bioclimatic variables mean at https://www.worldclim.org/data/bioclim.html.)

#element 0 is the intercept (predicted value if all climate variables are zero)
summary(m2)
## 
## Call:
## glm(formula = pb ~ bio1 + bio2 + bio3 + bio4 + bio5 + bio6 + 
##     bio7 + bio8 + bio9 + bio10 + bio11 + bio12 + bio13 + bio14 + 
##     bio15 + bio16 + bio17 + bio18 + bio19, data = sdmdata)
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.360e+00  4.521e-01   7.432 1.77e-13 ***
## bio1        -2.527e-02  9.514e-02  -0.266 0.790550    
## bio2         1.118e-01  5.008e-02   2.231 0.025792 *  
## bio3        -3.594e-02  6.001e-03  -5.989 2.62e-09 ***
## bio4        -2.665e-03  3.699e-03  -0.721 0.471207    
## bio5        -5.281e+04  3.093e+04  -1.707 0.087973 .  
## bio6         5.281e+04  3.093e+04   1.707 0.087973 .  
## bio7         5.281e+04  3.093e+04   1.707 0.087974 .  
## bio8         5.273e-02  1.907e-02   2.765 0.005758 ** 
## bio9         3.401e-02  2.212e-02   1.538 0.124339    
## bio10        4.933e-02  1.471e-01   0.335 0.737404    
## bio11       -1.398e-01  1.696e-01  -0.824 0.409858    
## bio12        3.228e-04  9.704e-05   3.326 0.000900 ***
## bio13        2.266e-03  6.382e-04   3.551 0.000395 ***
## bio14        3.021e-03  1.711e-03   1.766 0.077550 .  
## bio15        3.303e-04  9.191e-04   0.359 0.719355    
## bio16       -9.717e-04  3.149e-04  -3.085 0.002069 ** 
## bio17       -1.414e-03  6.490e-04  -2.179 0.029487 *  
## bio18        1.759e-04  1.158e-04   1.520 0.128811    
## bio19        1.381e-04  6.330e-05   2.182 0.029245 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 0.1768623)
## 
##     Null deviance: 358.97  on 1559  degrees of freedom
## Residual deviance: 272.37  on 1540  degrees of freedom
##   (1 observation deleted due to missingness)
## AIC: 1746.4
## 
## Number of Fisher Scoring iterations: 2
pvalues <- summary(m2)$coefficients[,4]
plot(pvalues, type = "n"); text(pvalues, labels = 0:19); segments(0, 0.05, 21, 0.05, col="red", lty=2)

# The red dashed line indicates p = 0.05, which is a traditional p-value 'cutoff' for assessing statistical  significance.

For my butterflies, the most predictive climate variables are:

  • bio2 (average daily temperature range) (somewhat less significant) *
  • bio3 (Isothermality (BIO2/BIO7) (×100)) ***
  • bio8 (Mean Temperature of Wettest Quarter) **
  • bio12 (Annual Precipitation) ***
  • bio13 (Precipitation of Wettest Month) ***
  • bio16 (Precipitation of Wettest Quarter) **
  • bio17 (Precipitation of Driest Quarter) *
  • bio19 (Precipitation of Coldest Quarter) *

But which direction does the influence of these variables go? Change the predictor names in “plotpredictors” to the most predictive climate variables for your species, then run the following code to generate box plots. Here, I’m just choosing the most significant variables: bio3, bio12, bio13. These will let you check whether higher or lower values of the predictor are associated with species occurrence.

plotpredictors <- c("bio3", "bio12", "bio13")
par(mfrow=c(1,length(plotpredictors)))
for (i in plotpredictors){
  boxplot(sdmdata[[i]] ~ sdmdata$pb, 
          xlab = "species presence", 
          ylab = "value of bioclim variable",
          main = i)
}

In these box plots, the points with x = 1 show the values of the predictor at the real occurrence points, while points with x = 0 show the values of the predictor variable at each of the fake absence points that we added. For my example species, compared to the absence points scattered across the rest of Central and South America, the presence points usually have a narrower range of values for these three significant variables. They prefer high isothermality (bio3, places with relatively high daily temperature fluctuations), high precipitation overall (bio12) and especially precipitation in the wettest month (bio13). In other words, it seems like blue morpho butterflies like living in the tropics because the tropics are wet, and they particularly like areas with fluctuating temperatures.


Question: Your species will certainly be different from mine. Which climate variables were most important to them? What do these variables mean?


Now to make the prediction map for the linear model:

q <- predict(climate,m2,progress = "text")
## |---------|---------|---------|---------|=========================================                                          

Save and reload it:

q <- reload_raster(q)

And plot it with the occurrence points:

plot(q, main = "Linear model prediction",
  xlim = range(df$lon,na.rm = T), 
  ylim = range(df$lat,na.rm = T))
points(df$lon,df$lat,col = rgb(1,0,0,0.5),pch = 20,cex = 0.5)

The difference between the bioclim model and the glm model, other than the scale of the suitability metric, is that the former doesn’t use the information about combinations of the different climate variables. You can see that the linear model did better at overlapping with the blue morpho’s real range – the areas with higher occurrence densities are more green.


Question: How did the linear model do with your species? Can you notice any differences between the Bioclim and linear models?


Comparing the different models

Higher value at a point on the map indicates that the model “thinks” the climate at that point is a good match for what the species likes (i.e. high habitat suitability). The maps we made show the locations of the actual occurrences, which lets you visually assess whether the real observations were found in places where the species is “expected” to live. Let’s try assessing this quantitatively now, by extracting the values of the model predictions at the observation locations and seeing whether each model is accurately predicting where the species has been observed.

modelmatch <- data.frame(bioclim = extract(p, df[,c("lon","lat")]),
                         linear = extract(q, df[,c("lon","lat")]))

Plot histograms of the results to compare them.

par(mfrow = c(2,1), mar = c(1,0,1,0))
hist(na.omit(modelmatch$bioclim/max(modelmatch$bioclim,na.rm = T)), col = "darkblue", main = "")
legend("topright",legend = names(modelmatch), col = c("darkblue","turquoise"), pch = 15)
hist(na.omit(modelmatch$linear/max(modelmatch$linear,na.rm = T)), col = "turquoise", main = "")

This plot summarizes the climate suitability values at the real occurrence points for each model. You can see that for Morpho menelaus, the values in the bioclim model are relatively lower, and the linear values are moderate. That is: the bioclim model did a worse job of predicting habitat suitability at the points where a blue morpho butterfly was seen, and the linear model did better (but maybe not great).


Question: Which modeling algorithm did the best job matching the true occurrence points for your species?


Save your plots, then save the data and unload:

save(raw_data, sfdata, df, m1, m2, p, q, sdmdata, modelmatch, file = here("climatemodel.RData"))
rm(raw_data, sfdata, df, m1, m2, p, q, sdmdata, modelmatch, climate, backgr, presvals, absvals)
gc()
##           used  (Mb) gc trigger   (Mb) limit (Mb)  max used   (Mb)
## Ncells 3026291 161.7   13351118  713.1         NA   8724651  466.0
## Vcells 6256975  47.8  476202612 3633.2      36864 744066580 5676.8

If you need to look at these again, run the following line:

load(here("climatemodel.RData"))

Part 2. Climate change and environmental variables

So far we’ve used only modern climate data and occurrences. However, we know that climate changes over time: 22,000 years ago we were in an ice age, and many species occurred in very different places than they are found today. Climate scientists have produced very detailed maps of the climatic conditions on Earth in the relatively recent past (especially for the last 3 million years or so). What can we do with estimates of past climate – that is, paleoclimate?

In this section you are going to be introduced to paleoclimate data sets. We’ll start by examining the mean annual temperature for your current location (UC Merced) in the past and present (and also for the projected future), and then you’ll do the same for another location that you choose yourself.

To load the paleoclimate and future climate data, run the following:

source("assemble-paleoclimate.R")

Examine the objects that it just loaded. These are four RasterStacks: lgm (Last Glacial Maximum, about 22,000 years ago), midH (Mid-Holocene, about 6000 years ago), modern (2000 CE), and future (projected climate for 2070 CE). Example code is below, for the LGM.

str(lgm)
## Formal class 'RasterStack' [package "raster"] with 12 slots
##   ..@ filename: chr ""
##   ..@ layers  :List of 19
##   .. ..$ :Formal class 'RasterLayer' [package "raster"] with 13 slots
##   .. .. .. ..@ file    :Formal class '.RasterFile' [package "raster"] with 13 slots
##   .. .. .. .. .. ..@ name        : chr "/Users/jblois/Documents/GitHub/biodata_shortcourse/development/paleoclimate/chelsa_LGM_v1_2B_r2_5m/bio_1.tif"
##   .. .. .. .. .. ..@ datanotation: chr "FLT4S"
##   .. .. .. .. .. ..@ byteorder   : chr "little"
##   .. .. .. .. .. ..@ nodatavalue : num -Inf
##   .. .. .. .. .. ..@ NAchanged   : logi FALSE
##   .. .. .. .. .. ..@ nbands      : int 1
##   .. .. .. .. .. ..@ bandorder   : chr "BIL"
##   .. .. .. .. .. ..@ offset      : int 0
##   .. .. .. .. .. ..@ toptobottom : logi TRUE
##   .. .. .. .. .. ..@ blockrows   : Named int 1
##   .. .. .. .. .. .. ..- attr(*, "names")= chr "rows"
##   .. .. .. .. .. ..@ blockcols   : Named int 8640
##   .. .. .. .. .. .. ..- attr(*, "names")= chr "cols"
##   .. .. .. .. .. ..@ driver      : chr "gdal"
##   .. .. .. .. .. ..@ open        : logi FALSE
##   .. .. .. ..@ data    :Formal class '.SingleLayerData' [package "raster"] with 13 slots
##   .. .. .. .. .. ..@ values    : logi(0) 
##   .. .. .. .. .. ..@ offset    : Named num 0
##   .. .. .. .. .. .. ..- attr(*, "names")= chr "offset"
##   .. .. .. .. .. ..@ gain      : Named num 1
##   .. .. .. .. .. .. ..- attr(*, "names")= chr "scale"
##   .. .. .. .. .. ..@ inmemory  : logi FALSE
##   .. .. .. .. .. ..@ fromdisk  : logi TRUE
##   .. .. .. .. .. ..@ isfactor  : logi FALSE
##   .. .. .. .. .. ..@ attributes: list()
##   .. .. .. .. .. ..@ haveminmax: logi TRUE
##   .. .. .. .. .. ..@ min       : num -80
##   .. .. .. .. .. ..@ max       : num 28.1
##   .. .. .. .. .. ..@ band      : int 1
##   .. .. .. .. .. ..@ unit      : chr ""
##   .. .. .. .. .. ..@ names     : chr "bio1"
##   .. .. .. ..@ legend  :Formal class '.RasterLegend' [package "raster"] with 5 slots
##   .. .. .. .. .. ..@ type      : chr(0) 
##   .. .. .. .. .. ..@ values    : logi(0) 
##   .. .. .. .. .. ..@ color     : logi(0) 
##   .. .. .. .. .. ..@ names     : logi(0) 
##   .. .. .. .. .. ..@ colortable: logi(0) 
##   .. .. .. ..@ title   : chr(0) 
##   .. .. .. ..@ extent  :Formal class 'Extent' [package "raster"] with 4 slots
##   .. .. .. .. .. ..@ xmin: num -180
##   .. .. .. .. .. ..@ xmax: num 180
##   .. .. .. .. .. ..@ ymin: num -90
##   .. .. .. .. .. ..@ ymax: num 90
##   .. .. .. ..@ rotated : logi FALSE
##   .. .. .. ..@ rotation:Formal class '.Rotation' [package "raster"] with 2 slots
##   .. .. .. .. .. ..@ geotrans: num(0) 
##   .. .. .. .. .. ..@ transfun:function ()  
##   .. .. .. ..@ ncols   : int 8640
##   .. .. .. ..@ nrows   : int 4320
##   .. .. .. ..@ crs     :Formal class 'CRS' [package "sp"] with 1 slot
##   .. .. .. .. .. ..@ projargs: chr NA
##   .. .. .. ..@ srs     : chr "+proj=longlat +datum=WGS84 +no_defs"
##   .. .. .. ..@ history : list()
##   .. .. .. ..@ z       : list()
##   .. ..$ :Formal class 'RasterLayer' [package "raster"] with 13 slots
##   .. .. .. ..@ file    :Formal class '.RasterFile' [package "raster"] with 13 slots
##   .. .. .. .. .. ..@ name        : chr "/Users/jblois/Documents/GitHub/biodata_shortcourse/development/paleoclimate/chelsa_LGM_v1_2B_r2_5m/bio_2.tif"
##   .. .. .. .. .. ..@ datanotation: chr "FLT4S"
##   .. .. .. .. .. ..@ byteorder   : chr "little"
##   .. .. .. .. .. ..@ nodatavalue : num -Inf
##   .. .. .. .. .. ..@ NAchanged   : logi FALSE
##   .. .. .. .. .. ..@ nbands      : int 1
##   .. .. .. .. .. ..@ bandorder   : chr "BIL"
##   .. .. .. .. .. ..@ offset      : int 0
##   .. .. .. .. .. ..@ toptobottom : logi TRUE
##   .. .. .. .. .. ..@ blockrows   : Named int 1
##   .. .. .. .. .. .. ..- attr(*, "names")= chr "rows"
##   .. .. .. .. .. ..@ blockcols   : Named int 8640
##   .. .. .. .. .. .. ..- attr(*, "names")= chr "cols"
##   .. .. .. .. .. ..@ driver      : chr "gdal"
##   .. .. .. .. .. ..@ open        : logi FALSE
##   .. .. .. ..@ data    :Formal class '.SingleLayerData' [package "raster"] with 13 slots
##   .. .. .. .. .. ..@ values    : logi(0) 
##   .. .. .. .. .. ..@ offset    : Named num 0
##   .. .. .. .. .. .. ..- attr(*, "names")= chr "offset"
##   .. .. .. .. .. ..@ gain      : Named num 1
##   .. .. .. .. .. .. ..- attr(*, "names")= chr "scale"
##   .. .. .. .. .. ..@ inmemory  : logi FALSE
##   .. .. .. .. .. ..@ fromdisk  : logi TRUE
##   .. .. .. .. .. ..@ isfactor  : logi FALSE
##   .. .. .. .. .. ..@ attributes: list()
##   .. .. .. .. .. ..@ haveminmax: logi TRUE
##   .. .. .. .. .. ..@ min       : num -1.4
##   .. .. .. .. .. ..@ max       : num 17.9
##   .. .. .. .. .. ..@ band      : int 1
##   .. .. .. .. .. ..@ unit      : chr ""
##   .. .. .. .. .. ..@ names     : chr "bio2"
##   .. .. .. ..@ legend  :Formal class '.RasterLegend' [package "raster"] with 5 slots
##   .. .. .. .. .. ..@ type      : chr(0) 
##   .. .. .. .. .. ..@ values    : logi(0) 
##   .. .. .. .. .. ..@ color     : logi(0) 
##   .. .. .. .. .. ..@ names     : logi(0) 
##   .. .. .. .. .. ..@ colortable: logi(0) 
##   .. .. .. ..@ title   : chr(0) 
##   .. .. .. ..@ extent  :Formal class 'Extent' [package "raster"] with 4 slots
##   .. .. .. .. .. ..@ xmin: num -180
##   .. .. .. .. .. ..@ xmax: num 180
##   .. .. .. .. .. ..@ ymin: num -90
##   .. .. .. .. .. ..@ ymax: num 90
##   .. .. .. ..@ rotated : logi FALSE
##   .. .. .. ..@ rotation:Formal class '.Rotation' [package "raster"] with 2 slots
##   .. .. .. .. .. ..@ geotrans: num(0) 
##   .. .. .. .. .. ..@ transfun:function ()  
##   .. .. .. ..@ ncols   : int 8640
##   .. .. .. ..@ nrows   : int 4320
##   .. .. .. ..@ crs     :Formal class 'CRS' [package "sp"] with 1 slot
##   .. .. .. .. .. ..@ projargs: chr NA
##   .. .. .. ..@ srs     : chr "+proj=longlat +datum=WGS84 +no_defs"
##   .. .. .. ..@ history : list()
##   .. .. .. ..@ z       : list()
##   .. ..$ :Formal class 'RasterLayer' [package "raster"] with 13 slots
##   .. .. .. ..@ file    :Formal class '.RasterFile' [package "raster"] with 13 slots
##   .. .. .. .. .. ..@ name        : chr "/Users/jblois/Documents/GitHub/biodata_shortcourse/development/paleoclimate/chelsa_LGM_v1_2B_r2_5m/bio_3.tif"
##   .. .. .. .. .. ..@ datanotation: chr "INT1U"
##   .. .. .. .. .. ..@ byteorder   : chr "little"
##   .. .. .. .. .. ..@ nodatavalue : num -Inf
##   .. .. .. .. .. ..@ NAchanged   : logi FALSE
##   .. .. .. .. .. ..@ nbands      : int 1
##   .. .. .. .. .. ..@ bandorder   : chr "BIL"
##   .. .. .. .. .. ..@ offset      : int 0
##   .. .. .. .. .. ..@ toptobottom : logi TRUE
##   .. .. .. .. .. ..@ blockrows   : Named int 1
##   .. .. .. .. .. .. ..- attr(*, "names")= chr "rows"
##   .. .. .. .. .. ..@ blockcols   : Named int 8640
##   .. .. .. .. .. .. ..- attr(*, "names")= chr "cols"
##   .. .. .. .. .. ..@ driver      : chr "gdal"
##   .. .. .. .. .. ..@ open        : logi FALSE
##   .. .. .. ..@ data    :Formal class '.SingleLayerData' [package "raster"] with 13 slots
##   .. .. .. .. .. ..@ values    : logi(0) 
##   .. .. .. .. .. ..@ offset    : Named num 0
##   .. .. .. .. .. .. ..- attr(*, "names")= chr "offset"
##   .. .. .. .. .. ..@ gain      : Named num 1
##   .. .. .. .. .. .. ..- attr(*, "names")= chr "scale"
##   .. .. .. .. .. ..@ inmemory  : logi FALSE
##   .. .. .. .. .. ..@ fromdisk  : logi TRUE
##   .. .. .. .. .. ..@ isfactor  : logi FALSE
##   .. .. .. .. .. ..@ attributes: list()
##   .. .. .. .. .. ..@ haveminmax: logi TRUE
##   .. .. .. .. .. ..@ min       : num 0
##   .. .. .. .. .. ..@ max       : num 255
##   .. .. .. .. .. ..@ band      : int 1
##   .. .. .. .. .. ..@ unit      : chr ""
##   .. .. .. .. .. ..@ names     : chr "bio3"
##   .. .. .. ..@ legend  :Formal class '.RasterLegend' [package "raster"] with 5 slots
##   .. .. .. .. .. ..@ type      : chr(0) 
##   .. .. .. .. .. ..@ values    : logi(0) 
##   .. .. .. .. .. ..@ color     : logi(0) 
##   .. .. .. .. .. ..@ names     : logi(0) 
##   .. .. .. .. .. ..@ colortable: logi(0) 
##   .. .. .. ..@ title   : chr(0) 
##   .. .. .. ..@ extent  :Formal class 'Extent' [package "raster"] with 4 slots
##   .. .. .. .. .. ..@ xmin: num -180
##   .. .. .. .. .. ..@ xmax: num 180
##   .. .. .. .. .. ..@ ymin: num -90
##   .. .. .. .. .. ..@ ymax: num 90
##   .. .. .. ..@ rotated : logi FALSE
##   .. .. .. ..@ rotation:Formal class '.Rotation' [package "raster"] with 2 slots
##   .. .. .. .. .. ..@ geotrans: num(0) 
##   .. .. .. .. .. ..@ transfun:function ()  
##   .. .. .. ..@ ncols   : int 8640
##   .. .. .. ..@ nrows   : int 4320
##   .. .. .. ..@ crs     :Formal class 'CRS' [package "sp"] with 1 slot
##   .. .. .. .. .. ..@ projargs: chr NA
##   .. .. .. ..@ srs     : chr "+proj=longlat +datum=WGS84 +no_defs"
##   .. .. .. ..@ history : list()
##   .. .. .. ..@ z       : list()
##   .. ..$ :Formal class 'RasterLayer' [package "raster"] with 13 slots
##   .. .. .. ..@ file    :Formal class '.RasterFile' [package "raster"] with 13 slots
##   .. .. .. .. .. ..@ name        : chr "/Users/jblois/Documents/GitHub/biodata_shortcourse/development/paleoclimate/chelsa_LGM_v1_2B_r2_5m/bio_4.tif"
##   .. .. .. .. .. ..@ datanotation: chr "FLT4S"
##   .. .. .. .. .. ..@ byteorder   : chr "little"
##   .. .. .. .. .. ..@ nodatavalue : num -Inf
##   .. .. .. .. .. ..@ NAchanged   : logi FALSE
##   .. .. .. .. .. ..@ nbands      : int 1
##   .. .. .. .. .. ..@ bandorder   : chr "BIL"
##   .. .. .. .. .. ..@ offset      : int 0
##   .. .. .. .. .. ..@ toptobottom : logi TRUE
##   .. .. .. .. .. ..@ blockrows   : Named int 1
##   .. .. .. .. .. .. ..- attr(*, "names")= chr "rows"
##   .. .. .. .. .. ..@ blockcols   : Named int 8640
##   .. .. .. .. .. .. ..- attr(*, "names")= chr "cols"
##   .. .. .. .. .. ..@ driver      : chr "gdal"
##   .. .. .. .. .. ..@ open        : logi FALSE
##   .. .. .. ..@ data    :Formal class '.SingleLayerData' [package "raster"] with 13 slots
##   .. .. .. .. .. ..@ values    : logi(0) 
##   .. .. .. .. .. ..@ offset    : Named num 0
##   .. .. .. .. .. .. ..- attr(*, "names")= chr "offset"
##   .. .. .. .. .. ..@ gain      : Named num 1
##   .. .. .. .. .. .. ..- attr(*, "names")= chr "scale"
##   .. .. .. .. .. ..@ inmemory  : logi FALSE
##   .. .. .. .. .. ..@ fromdisk  : logi TRUE
##   .. .. .. .. .. ..@ isfactor  : logi FALSE
##   .. .. .. .. .. ..@ attributes: list()
##   .. .. .. .. .. ..@ haveminmax: logi TRUE
##   .. .. .. .. .. ..@ min       : num 6.1
##   .. .. .. .. .. ..@ max       : num 2200
##   .. .. .. .. .. ..@ band      : int 1
##   .. .. .. .. .. ..@ unit      : chr ""
##   .. .. .. .. .. ..@ names     : chr "bio4"
##   .. .. .. ..@ legend  :Formal class '.RasterLegend' [package "raster"] with 5 slots
##   .. .. .. .. .. ..@ type      : chr(0) 
##   .. .. .. .. .. ..@ values    : logi(0) 
##   .. .. .. .. .. ..@ color     : logi(0) 
##   .. .. .. .. .. ..@ names     : logi(0) 
##   .. .. .. .. .. ..@ colortable: logi(0) 
##   .. .. .. ..@ title   : chr(0) 
##   .. .. .. ..@ extent  :Formal class 'Extent' [package "raster"] with 4 slots
##   .. .. .. .. .. ..@ xmin: num -180
##   .. .. .. .. .. ..@ xmax: num 180
##   .. .. .. .. .. ..@ ymin: num -90
##   .. .. .. .. .. ..@ ymax: num 90
##   .. .. .. ..@ rotated : logi FALSE
##   .. .. .. ..@ rotation:Formal class '.Rotation' [package "raster"] with 2 slots
##   .. .. .. .. .. ..@ geotrans: num(0) 
##   .. .. .. .. .. ..@ transfun:function ()  
##   .. .. .. ..@ ncols   : int 8640
##   .. .. .. ..@ nrows   : int 4320
##   .. .. .. ..@ crs     :Formal class 'CRS' [package "sp"] with 1 slot
##   .. .. .. .. .. ..@ projargs: chr NA
##   .. .. .. ..@ srs     : chr "+proj=longlat +datum=WGS84 +no_defs"
##   .. .. .. ..@ history : list()
##   .. .. .. ..@ z       : list()
##   .. ..$ :Formal class 'RasterLayer' [package "raster"] with 13 slots
##   .. .. .. ..@ file    :Formal class '.RasterFile' [package "raster"] with 13 slots
##   .. .. .. .. .. ..@ name        : chr "/Users/jblois/Documents/GitHub/biodata_shortcourse/development/paleoclimate/chelsa_LGM_v1_2B_r2_5m/bio_5.tif"
##   .. .. .. .. .. ..@ datanotation: chr "FLT4S"
##   .. .. .. .. .. ..@ byteorder   : chr "little"
##   .. .. .. .. .. ..@ nodatavalue : num -Inf
##   .. .. .. .. .. ..@ NAchanged   : logi FALSE
##   .. .. .. .. .. ..@ nbands      : int 1
##   .. .. .. .. .. ..@ bandorder   : chr "BIL"
##   .. .. .. .. .. ..@ offset      : int 0
##   .. .. .. .. .. ..@ toptobottom : logi TRUE
##   .. .. .. .. .. ..@ blockrows   : Named int 1
##   .. .. .. .. .. .. ..- attr(*, "names")= chr "rows"
##   .. .. .. .. .. ..@ blockcols   : Named int 8640
##   .. .. .. .. .. .. ..- attr(*, "names")= chr "cols"
##   .. .. .. .. .. ..@ driver      : chr "gdal"
##   .. .. .. .. .. ..@ open        : logi FALSE
##   .. .. .. ..@ data    :Formal class '.SingleLayerData' [package "raster"] with 13 slots
##   .. .. .. .. .. ..@ values    : logi(0) 
##   .. .. .. .. .. ..@ offset    : Named num 0
##   .. .. .. .. .. .. ..- attr(*, "names")= chr "offset"
##   .. .. .. .. .. ..@ gain      : Named num 1
##   .. .. .. .. .. .. ..- attr(*, "names")= chr "scale"
##   .. .. .. .. .. ..@ inmemory  : logi FALSE
##   .. .. .. .. .. ..@ fromdisk  : logi TRUE
##   .. .. .. .. .. ..@ isfactor  : logi FALSE
##   .. .. .. .. .. ..@ attributes: list()
##   .. .. .. .. .. ..@ haveminmax: logi TRUE
##   .. .. .. .. .. ..@ min       : num -55.3
##   .. .. .. .. .. ..@ max       : num 44.2
##   .. .. .. .. .. ..@ band      : int 1
##   .. .. .. .. .. ..@ unit      : chr ""
##   .. .. .. .. .. ..@ names     : chr "bio5"
##   .. .. .. ..@ legend  :Formal class '.RasterLegend' [package "raster"] with 5 slots
##   .. .. .. .. .. ..@ type      : chr(0) 
##   .. .. .. .. .. ..@ values    : logi(0) 
##   .. .. .. .. .. ..@ color     : logi(0) 
##   .. .. .. .. .. ..@ names     : logi(0) 
##   .. .. .. .. .. ..@ colortable: logi(0) 
##   .. .. .. ..@ title   : chr(0) 
##   .. .. .. ..@ extent  :Formal class 'Extent' [package "raster"] with 4 slots
##   .. .. .. .. .. ..@ xmin: num -180
##   .. .. .. .. .. ..@ xmax: num 180
##   .. .. .. .. .. ..@ ymin: num -90
##   .. .. .. .. .. ..@ ymax: num 90
##   .. .. .. ..@ rotated : logi FALSE
##   .. .. .. ..@ rotation:Formal class '.Rotation' [package "raster"] with 2 slots
##   .. .. .. .. .. ..@ geotrans: num(0) 
##   .. .. .. .. .. ..@ transfun:function ()  
##   .. .. .. ..@ ncols   : int 8640
##   .. .. .. ..@ nrows   : int 4320
##   .. .. .. ..@ crs     :Formal class 'CRS' [package "sp"] with 1 slot
##   .. .. .. .. .. ..@ projargs: chr NA
##   .. .. .. ..@ srs     : chr "+proj=longlat +datum=WGS84 +no_defs"
##   .. .. .. ..@ history : list()
##   .. .. .. ..@ z       : list()
##   .. ..$ :Formal class 'RasterLayer' [package "raster"] with 13 slots
##   .. .. .. ..@ file    :Formal class '.RasterFile' [package "raster"] with 13 slots
##   .. .. .. .. .. ..@ name        : chr "/Users/jblois/Documents/GitHub/biodata_shortcourse/development/paleoclimate/chelsa_LGM_v1_2B_r2_5m/bio_6.tif"
##   .. .. .. .. .. ..@ datanotation: chr "INT2S"
##   .. .. .. .. .. ..@ byteorder   : chr "little"
##   .. .. .. .. .. ..@ nodatavalue : num -Inf
##   .. .. .. .. .. ..@ NAchanged   : logi FALSE
##   .. .. .. .. .. ..@ nbands      : int 1
##   .. .. .. .. .. ..@ bandorder   : chr "BIL"
##   .. .. .. .. .. ..@ offset      : int 0
##   .. .. .. .. .. ..@ toptobottom : logi TRUE
##   .. .. .. .. .. ..@ blockrows   : Named int 1
##   .. .. .. .. .. .. ..- attr(*, "names")= chr "rows"
##   .. .. .. .. .. ..@ blockcols   : Named int 8640
##   .. .. .. .. .. .. ..- attr(*, "names")= chr "cols"
##   .. .. .. .. .. ..@ driver      : chr "gdal"
##   .. .. .. .. .. ..@ open        : logi FALSE
##   .. .. .. ..@ data    :Formal class '.SingleLayerData' [package "raster"] with 13 slots
##   .. .. .. .. .. ..@ values    : logi(0) 
##   .. .. .. .. .. ..@ offset    : Named num 0
##   .. .. .. .. .. .. ..- attr(*, "names")= chr "offset"
##   .. .. .. .. .. ..@ gain      : Named num 1
##   .. .. .. .. .. .. ..- attr(*, "names")= chr "scale"
##   .. .. .. .. .. ..@ inmemory  : logi FALSE
##   .. .. .. .. .. ..@ fromdisk  : logi TRUE
##   .. .. .. .. .. ..@ isfactor  : logi FALSE
##   .. .. .. .. .. ..@ attributes: list()
##   .. .. .. .. .. ..@ haveminmax: logi TRUE
##   .. .. .. .. .. ..@ min       : num -991
##   .. .. .. .. .. ..@ max       : num 250
##   .. .. .. .. .. ..@ band      : int 1
##   .. .. .. .. .. ..@ unit      : chr ""
##   .. .. .. .. .. ..@ names     : chr "bio6"
##   .. .. .. ..@ legend  :Formal class '.RasterLegend' [package "raster"] with 5 slots
##   .. .. .. .. .. ..@ type      : chr(0) 
##   .. .. .. .. .. ..@ values    : logi(0) 
##   .. .. .. .. .. ..@ color     : logi(0) 
##   .. .. .. .. .. ..@ names     : logi(0) 
##   .. .. .. .. .. ..@ colortable: logi(0) 
##   .. .. .. ..@ title   : chr(0) 
##   .. .. .. ..@ extent  :Formal class 'Extent' [package "raster"] with 4 slots
##   .. .. .. .. .. ..@ xmin: num -180
##   .. .. .. .. .. ..@ xmax: num 180
##   .. .. .. .. .. ..@ ymin: num -90
##   .. .. .. .. .. ..@ ymax: num 90
##   .. .. .. ..@ rotated : logi FALSE
##   .. .. .. ..@ rotation:Formal class '.Rotation' [package "raster"] with 2 slots
##   .. .. .. .. .. ..@ geotrans: num(0) 
##   .. .. .. .. .. ..@ transfun:function ()  
##   .. .. .. ..@ ncols   : int 8640
##   .. .. .. ..@ nrows   : int 4320
##   .. .. .. ..@ crs     :Formal class 'CRS' [package "sp"] with 1 slot
##   .. .. .. .. .. ..@ projargs: chr NA
##   .. .. .. ..@ srs     : chr "+proj=longlat +datum=WGS84 +no_defs"
##   .. .. .. ..@ history : list()
##   .. .. .. ..@ z       : list()
##   .. ..$ :Formal class 'RasterLayer' [package "raster"] with 13 slots
##   .. .. .. ..@ file    :Formal class '.RasterFile' [package "raster"] with 13 slots
##   .. .. .. .. .. ..@ name        : chr "/Users/jblois/Documents/GitHub/biodata_shortcourse/development/paleoclimate/chelsa_LGM_v1_2B_r2_5m/bio_7.tif"
##   .. .. .. .. .. ..@ datanotation: chr "FLT4S"
##   .. .. .. .. .. ..@ byteorder   : chr "little"
##   .. .. .. .. .. ..@ nodatavalue : num -Inf
##   .. .. .. .. .. ..@ NAchanged   : logi FALSE
##   .. .. .. .. .. ..@ nbands      : int 1
##   .. .. .. .. .. ..@ bandorder   : chr "BIL"
##   .. .. .. .. .. ..@ offset      : int 0
##   .. .. .. .. .. ..@ toptobottom : logi TRUE
##   .. .. .. .. .. ..@ blockrows   : Named int 1
##   .. .. .. .. .. .. ..- attr(*, "names")= chr "rows"
##   .. .. .. .. .. ..@ blockcols   : Named int 8640
##   .. .. .. .. .. .. ..- attr(*, "names")= chr "cols"
##   .. .. .. .. .. ..@ driver      : chr "gdal"
##   .. .. .. .. .. ..@ open        : logi FALSE
##   .. .. .. ..@ data    :Formal class '.SingleLayerData' [package "raster"] with 13 slots
##   .. .. .. .. .. ..@ values    : logi(0) 
##   .. .. .. .. .. ..@ offset    : Named num 0
##   .. .. .. .. .. .. ..- attr(*, "names")= chr "offset"
##   .. .. .. .. .. ..@ gain      : Named num 1
##   .. .. .. .. .. .. ..- attr(*, "names")= chr "scale"
##   .. .. .. .. .. ..@ inmemory  : logi FALSE
##   .. .. .. .. .. ..@ fromdisk  : logi TRUE
##   .. .. .. .. .. ..@ isfactor  : logi FALSE
##   .. .. .. .. .. ..@ attributes: list()
##   .. .. .. .. .. ..@ haveminmax: logi TRUE
##   .. .. .. .. .. ..@ min       : num 0.8
##   .. .. .. .. .. ..@ max       : num 73.4
##   .. .. .. .. .. ..@ band      : int 1
##   .. .. .. .. .. ..@ unit      : chr ""
##   .. .. .. .. .. ..@ names     : chr "bio7"
##   .. .. .. ..@ legend  :Formal class '.RasterLegend' [package "raster"] with 5 slots
##   .. .. .. .. .. ..@ type      : chr(0) 
##   .. .. .. .. .. ..@ values    : logi(0) 
##   .. .. .. .. .. ..@ color     : logi(0) 
##   .. .. .. .. .. ..@ names     : logi(0) 
##   .. .. .. .. .. ..@ colortable: logi(0) 
##   .. .. .. ..@ title   : chr(0) 
##   .. .. .. ..@ extent  :Formal class 'Extent' [package "raster"] with 4 slots
##   .. .. .. .. .. ..@ xmin: num -180
##   .. .. .. .. .. ..@ xmax: num 180
##   .. .. .. .. .. ..@ ymin: num -90
##   .. .. .. .. .. ..@ ymax: num 90
##   .. .. .. ..@ rotated : logi FALSE
##   .. .. .. ..@ rotation:Formal class '.Rotation' [package "raster"] with 2 slots
##   .. .. .. .. .. ..@ geotrans: num(0) 
##   .. .. .. .. .. ..@ transfun:function ()  
##   .. .. .. ..@ ncols   : int 8640
##   .. .. .. ..@ nrows   : int 4320
##   .. .. .. ..@ crs     :Formal class 'CRS' [package "sp"] with 1 slot
##   .. .. .. .. .. ..@ projargs: chr NA
##   .. .. .. ..@ srs     : chr "+proj=longlat +datum=WGS84 +no_defs"
##   .. .. .. ..@ history : list()
##   .. .. .. ..@ z       : list()
##   .. ..$ :Formal class 'RasterLayer' [package "raster"] with 13 slots
##   .. .. .. ..@ file    :Formal class '.RasterFile' [package "raster"] with 13 slots
##   .. .. .. .. .. ..@ name        : chr "/Users/jblois/Documents/GitHub/biodata_shortcourse/development/paleoclimate/chelsa_LGM_v1_2B_r2_5m/bio_8.tif"
##   .. .. .. .. .. ..@ datanotation: chr "FLT4S"
##   .. .. .. .. .. ..@ byteorder   : chr "little"
##   .. .. .. .. .. ..@ nodatavalue : num -Inf
##   .. .. .. .. .. ..@ NAchanged   : logi FALSE
##   .. .. .. .. .. ..@ nbands      : int 1
##   .. .. .. .. .. ..@ bandorder   : chr "BIL"
##   .. .. .. .. .. ..@ offset      : int 0
##   .. .. .. .. .. ..@ toptobottom : logi TRUE
##   .. .. .. .. .. ..@ blockrows   : Named int 1
##   .. .. .. .. .. .. ..- attr(*, "names")= chr "rows"
##   .. .. .. .. .. ..@ blockcols   : Named int 8640
##   .. .. .. .. .. .. ..- attr(*, "names")= chr "cols"
##   .. .. .. .. .. ..@ driver      : chr "gdal"
##   .. .. .. .. .. ..@ open        : logi FALSE
##   .. .. .. ..@ data    :Formal class '.SingleLayerData' [package "raster"] with 13 slots
##   .. .. .. .. .. ..@ values    : logi(0) 
##   .. .. .. .. .. ..@ offset    : Named num 0
##   .. .. .. .. .. .. ..- attr(*, "names")= chr "offset"
##   .. .. .. .. .. ..@ gain      : Named num 1
##   .. .. .. .. .. .. ..- attr(*, "names")= chr "scale"
##   .. .. .. .. .. ..@ inmemory  : logi FALSE
##   .. .. .. .. .. ..@ fromdisk  : logi TRUE
##   .. .. .. .. .. ..@ isfactor  : logi FALSE
##   .. .. .. .. .. ..@ attributes: list()
##   .. .. .. .. .. ..@ haveminmax: logi TRUE
##   .. .. .. .. .. ..@ min       : num -87.7
##   .. .. .. .. .. ..@ max       : num 33.3
##   .. .. .. .. .. ..@ band      : int 1
##   .. .. .. .. .. ..@ unit      : chr ""
##   .. .. .. .. .. ..@ names     : chr "bio8"
##   .. .. .. ..@ legend  :Formal class '.RasterLegend' [package "raster"] with 5 slots
##   .. .. .. .. .. ..@ type      : chr(0) 
##   .. .. .. .. .. ..@ values    : logi(0) 
##   .. .. .. .. .. ..@ color     : logi(0) 
##   .. .. .. .. .. ..@ names     : logi(0) 
##   .. .. .. .. .. ..@ colortable: logi(0) 
##   .. .. .. ..@ title   : chr(0) 
##   .. .. .. ..@ extent  :Formal class 'Extent' [package "raster"] with 4 slots
##   .. .. .. .. .. ..@ xmin: num -180
##   .. .. .. .. .. ..@ xmax: num 180
##   .. .. .. .. .. ..@ ymin: num -90
##   .. .. .. .. .. ..@ ymax: num 90
##   .. .. .. ..@ rotated : logi FALSE
##   .. .. .. ..@ rotation:Formal class '.Rotation' [package "raster"] with 2 slots
##   .. .. .. .. .. ..@ geotrans: num(0) 
##   .. .. .. .. .. ..@ transfun:function ()  
##   .. .. .. ..@ ncols   : int 8640
##   .. .. .. ..@ nrows   : int 4320
##   .. .. .. ..@ crs     :Formal class 'CRS' [package "sp"] with 1 slot
##   .. .. .. .. .. ..@ projargs: chr NA
##   .. .. .. ..@ srs     : chr "+proj=longlat +datum=WGS84 +no_defs"
##   .. .. .. ..@ history : list()
##   .. .. .. ..@ z       : list()
##   .. ..$ :Formal class 'RasterLayer' [package "raster"] with 13 slots
##   .. .. .. ..@ file    :Formal class '.RasterFile' [package "raster"] with 13 slots
##   .. .. .. .. .. ..@ name        : chr "/Users/jblois/Documents/GitHub/biodata_shortcourse/development/paleoclimate/chelsa_LGM_v1_2B_r2_5m/bio_9.tif"
##   .. .. .. .. .. ..@ datanotation: chr "FLT4S"
##   .. .. .. .. .. ..@ byteorder   : chr "little"
##   .. .. .. .. .. ..@ nodatavalue : num -Inf
##   .. .. .. .. .. ..@ NAchanged   : logi FALSE
##   .. .. .. .. .. ..@ nbands      : int 1
##   .. .. .. .. .. ..@ bandorder   : chr "BIL"
##   .. .. .. .. .. ..@ offset      : int 0
##   .. .. .. .. .. ..@ toptobottom : logi TRUE
##   .. .. .. .. .. ..@ blockrows   : Named int 1
##   .. .. .. .. .. .. ..- attr(*, "names")= chr "rows"
##   .. .. .. .. .. ..@ blockcols   : Named int 8640
##   .. .. .. .. .. .. ..- attr(*, "names")= chr "cols"
##   .. .. .. .. .. ..@ driver      : chr "gdal"
##   .. .. .. .. .. ..@ open        : logi FALSE
##   .. .. .. ..@ data    :Formal class '.SingleLayerData' [package "raster"] with 13 slots
##   .. .. .. .. .. ..@ values    : logi(0) 
##   .. .. .. .. .. ..@ offset    : Named num 0
##   .. .. .. .. .. .. ..- attr(*, "names")= chr "offset"
##   .. .. .. .. .. ..@ gain      : Named num 1
##   .. .. .. .. .. .. ..- attr(*, "names")= chr "scale"
##   .. .. .. .. .. ..@ inmemory  : logi FALSE
##   .. .. .. .. .. ..@ fromdisk  : logi TRUE
##   .. .. .. .. .. ..@ isfactor  : logi FALSE
##   .. .. .. .. .. ..@ attributes: list()
##   .. .. .. .. .. ..@ haveminmax: logi TRUE
##   .. .. .. .. .. ..@ min       : num -85.7
##   .. .. .. .. .. ..@ max       : num 36.2
##   .. .. .. .. .. ..@ band      : int 1
##   .. .. .. .. .. ..@ unit      : chr ""
##   .. .. .. .. .. ..@ names     : chr "bio9"
##   .. .. .. ..@ legend  :Formal class '.RasterLegend' [package "raster"] with 5 slots
##   .. .. .. .. .. ..@ type      : chr(0) 
##   .. .. .. .. .. ..@ values    : logi(0) 
##   .. .. .. .. .. ..@ color     : logi(0) 
##   .. .. .. .. .. ..@ names     : logi(0) 
##   .. .. .. .. .. ..@ colortable: logi(0) 
##   .. .. .. ..@ title   : chr(0) 
##   .. .. .. ..@ extent  :Formal class 'Extent' [package "raster"] with 4 slots
##   .. .. .. .. .. ..@ xmin: num -180
##   .. .. .. .. .. ..@ xmax: num 180
##   .. .. .. .. .. ..@ ymin: num -90
##   .. .. .. .. .. ..@ ymax: num 90
##   .. .. .. ..@ rotated : logi FALSE
##   .. .. .. ..@ rotation:Formal class '.Rotation' [package "raster"] with 2 slots
##   .. .. .. .. .. ..@ geotrans: num(0) 
##   .. .. .. .. .. ..@ transfun:function ()  
##   .. .. .. ..@ ncols   : int 8640
##   .. .. .. ..@ nrows   : int 4320
##   .. .. .. ..@ crs     :Formal class 'CRS' [package "sp"] with 1 slot
##   .. .. .. .. .. ..@ projargs: chr NA
##   .. .. .. ..@ srs     : chr "+proj=longlat +datum=WGS84 +no_defs"
##   .. .. .. ..@ history : list()
##   .. .. .. ..@ z       : list()
##   .. ..$ :Formal class 'RasterLayer' [package "raster"] with 13 slots
##   .. .. .. ..@ file    :Formal class '.RasterFile' [package "raster"] with 13 slots
##   .. .. .. .. .. ..@ name        : chr "/Users/jblois/Documents/GitHub/biodata_shortcourse/development/paleoclimate/chelsa_LGM_v1_2B_r2_5m/bio_10.tif"
##   .. .. .. .. .. ..@ datanotation: chr "FLT4S"
##   .. .. .. .. .. ..@ byteorder   : chr "little"
##   .. .. .. .. .. ..@ nodatavalue : num -Inf
##   .. .. .. .. .. ..@ NAchanged   : logi FALSE
##   .. .. .. .. .. ..@ nbands      : int 1
##   .. .. .. .. .. ..@ bandorder   : chr "BIL"
##   .. .. .. .. .. ..@ offset      : int 0
##   .. .. .. .. .. ..@ toptobottom : logi TRUE
##   .. .. .. .. .. ..@ blockrows   : Named int 1
##   .. .. .. .. .. .. ..- attr(*, "names")= chr "rows"
##   .. .. .. .. .. ..@ blockcols   : Named int 8640
##   .. .. .. .. .. .. ..- attr(*, "names")= chr "cols"
##   .. .. .. .. .. ..@ driver      : chr "gdal"
##   .. .. .. .. .. ..@ open        : logi FALSE
##   .. .. .. ..@ data    :Formal class '.SingleLayerData' [package "raster"] with 13 slots
##   .. .. .. .. .. ..@ values    : logi(0) 
##   .. .. .. .. .. ..@ offset    : Named num 0
##   .. .. .. .. .. .. ..- attr(*, "names")= chr "offset"
##   .. .. .. .. .. ..@ gain      : Named num 1
##   .. .. .. .. .. .. ..- attr(*, "names")= chr "scale"
##   .. .. .. .. .. ..@ inmemory  : logi FALSE
##   .. .. .. .. .. ..@ fromdisk  : logi TRUE
##   .. .. .. .. .. ..@ isfactor  : logi FALSE
##   .. .. .. .. .. ..@ attributes: list()
##   .. .. .. .. .. ..@ haveminmax: logi TRUE
##   .. .. .. .. .. ..@ min       : num -57.8
##   .. .. .. .. .. ..@ max       : num 36.2
##   .. .. .. .. .. ..@ band      : int 1
##   .. .. .. .. .. ..@ unit      : chr ""
##   .. .. .. .. .. ..@ names     : chr "bio10"
##   .. .. .. ..@ legend  :Formal class '.RasterLegend' [package "raster"] with 5 slots
##   .. .. .. .. .. ..@ type      : chr(0) 
##   .. .. .. .. .. ..@ values    : logi(0) 
##   .. .. .. .. .. ..@ color     : logi(0) 
##   .. .. .. .. .. ..@ names     : logi(0) 
##   .. .. .. .. .. ..@ colortable: logi(0) 
##   .. .. .. ..@ title   : chr(0) 
##   .. .. .. ..@ extent  :Formal class 'Extent' [package "raster"] with 4 slots
##   .. .. .. .. .. ..@ xmin: num -180
##   .. .. .. .. .. ..@ xmax: num 180
##   .. .. .. .. .. ..@ ymin: num -90
##   .. .. .. .. .. ..@ ymax: num 90
##   .. .. .. ..@ rotated : logi FALSE
##   .. .. .. ..@ rotation:Formal class '.Rotation' [package "raster"] with 2 slots
##   .. .. .. .. .. ..@ geotrans: num(0) 
##   .. .. .. .. .. ..@ transfun:function ()  
##   .. .. .. ..@ ncols   : int 8640
##   .. .. .. ..@ nrows   : int 4320
##   .. .. .. ..@ crs     :Formal class 'CRS' [package "sp"] with 1 slot
##   .. .. .. .. .. ..@ projargs: chr NA
##   .. .. .. ..@ srs     : chr "+proj=longlat +datum=WGS84 +no_defs"
##   .. .. .. ..@ history : list()
##   .. .. .. ..@ z       : list()
##   .. ..$ :Formal class 'RasterLayer' [package "raster"] with 13 slots
##   .. .. .. ..@ file    :Formal class '.RasterFile' [package "raster"] with 13 slots
##   .. .. .. .. .. ..@ name        : chr "/Users/jblois/Documents/GitHub/biodata_shortcourse/development/paleoclimate/chelsa_LGM_v1_2B_r2_5m/bio_11.tif"
##   .. .. .. .. .. ..@ datanotation: chr "FLT4S"
##   .. .. .. .. .. ..@ byteorder   : chr "little"
##   .. .. .. .. .. ..@ nodatavalue : num -Inf
##   .. .. .. .. .. ..@ NAchanged   : logi FALSE
##   .. .. .. .. .. ..@ nbands      : int 1
##   .. .. .. .. .. ..@ bandorder   : chr "BIL"
##   .. .. .. .. .. ..@ offset      : int 0
##   .. .. .. .. .. ..@ toptobottom : logi TRUE
##   .. .. .. .. .. ..@ blockrows   : Named int 1
##   .. .. .. .. .. .. ..- attr(*, "names")= chr "rows"
##   .. .. .. .. .. ..@ blockcols   : Named int 8640
##   .. .. .. .. .. .. ..- attr(*, "names")= chr "cols"
##   .. .. .. .. .. ..@ driver      : chr "gdal"
##   .. .. .. .. .. ..@ open        : logi FALSE
##   .. .. .. ..@ data    :Formal class '.SingleLayerData' [package "raster"] with 13 slots
##   .. .. .. .. .. ..@ values    : logi(0) 
##   .. .. .. .. .. ..@ offset    : Named num 0
##   .. .. .. .. .. .. ..- attr(*, "names")= chr "offset"
##   .. .. .. .. .. ..@ gain      : Named num 1
##   .. .. .. .. .. .. ..- attr(*, "names")= chr "scale"
##   .. .. .. .. .. ..@ inmemory  : logi FALSE
##   .. .. .. .. .. ..@ fromdisk  : logi TRUE
##   .. .. .. .. .. ..@ isfactor  : logi FALSE
##   .. .. .. .. .. ..@ attributes: list()
##   .. .. .. .. .. ..@ haveminmax: logi TRUE
##   .. .. .. .. .. ..@ min       : num -97.2
##   .. .. .. .. .. ..@ max       : num 27.1
##   .. .. .. .. .. ..@ band      : int 1
##   .. .. .. .. .. ..@ unit      : chr ""
##   .. .. .. .. .. ..@ names     : chr "bio11"
##   .. .. .. ..@ legend  :Formal class '.RasterLegend' [package "raster"] with 5 slots
##   .. .. .. .. .. ..@ type      : chr(0) 
##   .. .. .. .. .. ..@ values    : logi(0) 
##   .. .. .. .. .. ..@ color     : logi(0) 
##   .. .. .. .. .. ..@ names     : logi(0) 
##   .. .. .. .. .. ..@ colortable: logi(0) 
##   .. .. .. ..@ title   : chr(0) 
##   .. .. .. ..@ extent  :Formal class 'Extent' [package "raster"] with 4 slots
##   .. .. .. .. .. ..@ xmin: num -180
##   .. .. .. .. .. ..@ xmax: num 180
##   .. .. .. .. .. ..@ ymin: num -90
##   .. .. .. .. .. ..@ ymax: num 90
##   .. .. .. ..@ rotated : logi FALSE
##   .. .. .. ..@ rotation:Formal class '.Rotation' [package "raster"] with 2 slots
##   .. .. .. .. .. ..@ geotrans: num(0) 
##   .. .. .. .. .. ..@ transfun:function ()  
##   .. .. .. ..@ ncols   : int 8640
##   .. .. .. ..@ nrows   : int 4320
##   .. .. .. ..@ crs     :Formal class 'CRS' [package "sp"] with 1 slot
##   .. .. .. .. .. ..@ projargs: chr NA
##   .. .. .. ..@ srs     : chr "+proj=longlat +datum=WGS84 +no_defs"
##   .. .. .. ..@ history : list()
##   .. .. .. ..@ z       : list()
##   .. ..$ :Formal class 'RasterLayer' [package "raster"] with 13 slots
##   .. .. .. ..@ file    :Formal class '.RasterFile' [package "raster"] with 13 slots
##   .. .. .. .. .. ..@ name        : chr "/Users/jblois/Documents/GitHub/biodata_shortcourse/development/paleoclimate/chelsa_LGM_v1_2B_r2_5m/bio_12.tif"
##   .. .. .. .. .. ..@ datanotation: chr "INT2U"
##   .. .. .. .. .. ..@ byteorder   : chr "little"
##   .. .. .. .. .. ..@ nodatavalue : num -Inf
##   .. .. .. .. .. ..@ NAchanged   : logi FALSE
##   .. .. .. .. .. ..@ nbands      : int 1
##   .. .. .. .. .. ..@ bandorder   : chr "BIL"
##   .. .. .. .. .. ..@ offset      : int 0
##   .. .. .. .. .. ..@ toptobottom : logi TRUE
##   .. .. .. .. .. ..@ blockrows   : Named int 1
##   .. .. .. .. .. .. ..- attr(*, "names")= chr "rows"
##   .. .. .. .. .. ..@ blockcols   : Named int 8640
##   .. .. .. .. .. .. ..- attr(*, "names")= chr "cols"
##   .. .. .. .. .. ..@ driver      : chr "gdal"
##   .. .. .. .. .. ..@ open        : logi FALSE
##   .. .. .. ..@ data    :Formal class '.SingleLayerData' [package "raster"] with 13 slots
##   .. .. .. .. .. ..@ values    : logi(0) 
##   .. .. .. .. .. ..@ offset    : Named num 0
##   .. .. .. .. .. .. ..- attr(*, "names")= chr "offset"
##   .. .. .. .. .. ..@ gain      : Named num 1
##   .. .. .. .. .. .. ..- attr(*, "names")= chr "scale"
##   .. .. .. .. .. ..@ inmemory  : logi FALSE
##   .. .. .. .. .. ..@ fromdisk  : logi TRUE
##   .. .. .. .. .. ..@ isfactor  : logi FALSE
##   .. .. .. .. .. ..@ attributes: list()
##   .. .. .. .. .. ..@ haveminmax: logi TRUE
##   .. .. .. .. .. ..@ min       : num 0
##   .. .. .. .. .. ..@ max       : num 3276
##   .. .. .. .. .. ..@ band      : int 1
##   .. .. .. .. .. ..@ unit      : chr ""
##   .. .. .. .. .. ..@ names     : chr "bio12"
##   .. .. .. ..@ legend  :Formal class '.RasterLegend' [package "raster"] with 5 slots
##   .. .. .. .. .. ..@ type      : chr(0) 
##   .. .. .. .. .. ..@ values    : logi(0) 
##   .. .. .. .. .. ..@ color     : logi(0) 
##   .. .. .. .. .. ..@ names     : logi(0) 
##   .. .. .. .. .. ..@ colortable: logi(0) 
##   .. .. .. ..@ title   : chr(0) 
##   .. .. .. ..@ extent  :Formal class 'Extent' [package "raster"] with 4 slots
##   .. .. .. .. .. ..@ xmin: num -180
##   .. .. .. .. .. ..@ xmax: num 180
##   .. .. .. .. .. ..@ ymin: num -90
##   .. .. .. .. .. ..@ ymax: num 90
##   .. .. .. ..@ rotated : logi FALSE
##   .. .. .. ..@ rotation:Formal class '.Rotation' [package "raster"] with 2 slots
##   .. .. .. .. .. ..@ geotrans: num(0) 
##   .. .. .. .. .. ..@ transfun:function ()  
##   .. .. .. ..@ ncols   : int 8640
##   .. .. .. ..@ nrows   : int 4320
##   .. .. .. ..@ crs     :Formal class 'CRS' [package "sp"] with 1 slot
##   .. .. .. .. .. ..@ projargs: chr NA
##   .. .. .. ..@ srs     : chr "+proj=longlat +datum=WGS84 +no_defs"
##   .. .. .. ..@ history : list()
##   .. .. .. ..@ z       : list()
##   .. ..$ :Formal class 'RasterLayer' [package "raster"] with 13 slots
##   .. .. .. ..@ file    :Formal class '.RasterFile' [package "raster"] with 13 slots
##   .. .. .. .. .. ..@ name        : chr "/Users/jblois/Documents/GitHub/biodata_shortcourse/development/paleoclimate/chelsa_LGM_v1_2B_r2_5m/bio_13.tif"
##   .. .. .. .. .. ..@ datanotation: chr "INT2U"
##   .. .. .. .. .. ..@ byteorder   : chr "little"
##   .. .. .. .. .. ..@ nodatavalue : num -Inf
##   .. .. .. .. .. ..@ NAchanged   : logi FALSE
##   .. .. .. .. .. ..@ nbands      : int 1
##   .. .. .. .. .. ..@ bandorder   : chr "BIL"
##   .. .. .. .. .. ..@ offset      : int 0
##   .. .. .. .. .. ..@ toptobottom : logi TRUE
##   .. .. .. .. .. ..@ blockrows   : Named int 1
##   .. .. .. .. .. .. ..- attr(*, "names")= chr "rows"
##   .. .. .. .. .. ..@ blockcols   : Named int 8640
##   .. .. .. .. .. .. ..- attr(*, "names")= chr "cols"
##   .. .. .. .. .. ..@ driver      : chr "gdal"
##   .. .. .. .. .. ..@ open        : logi FALSE
##   .. .. .. ..@ data    :Formal class '.SingleLayerData' [package "raster"] with 13 slots
##   .. .. .. .. .. ..@ values    : logi(0) 
##   .. .. .. .. .. ..@ offset    : Named num 0
##   .. .. .. .. .. .. ..- attr(*, "names")= chr "offset"
##   .. .. .. .. .. ..@ gain      : Named num 1
##   .. .. .. .. .. .. ..- attr(*, "names")= chr "scale"
##   .. .. .. .. .. ..@ inmemory  : logi FALSE
##   .. .. .. .. .. ..@ fromdisk  : logi TRUE
##   .. .. .. .. .. ..@ isfactor  : logi FALSE
##   .. .. .. .. .. ..@ attributes: list()
##   .. .. .. .. .. ..@ haveminmax: logi TRUE
##   .. .. .. .. .. ..@ min       : num 0
##   .. .. .. .. .. ..@ max       : num 2570
##   .. .. .. .. .. ..@ band      : int 1
##   .. .. .. .. .. ..@ unit      : chr ""
##   .. .. .. .. .. ..@ names     : chr "bio13"
##   .. .. .. ..@ legend  :Formal class '.RasterLegend' [package "raster"] with 5 slots
##   .. .. .. .. .. ..@ type      : chr(0) 
##   .. .. .. .. .. ..@ values    : logi(0) 
##   .. .. .. .. .. ..@ color     : logi(0) 
##   .. .. .. .. .. ..@ names     : logi(0) 
##   .. .. .. .. .. ..@ colortable: logi(0) 
##   .. .. .. ..@ title   : chr(0) 
##   .. .. .. ..@ extent  :Formal class 'Extent' [package "raster"] with 4 slots
##   .. .. .. .. .. ..@ xmin: num -180
##   .. .. .. .. .. ..@ xmax: num 180
##   .. .. .. .. .. ..@ ymin: num -90
##   .. .. .. .. .. ..@ ymax: num 90
##   .. .. .. ..@ rotated : logi FALSE
##   .. .. .. ..@ rotation:Formal class '.Rotation' [package "raster"] with 2 slots
##   .. .. .. .. .. ..@ geotrans: num(0) 
##   .. .. .. .. .. ..@ transfun:function ()  
##   .. .. .. ..@ ncols   : int 8640
##   .. .. .. ..@ nrows   : int 4320
##   .. .. .. ..@ crs     :Formal class 'CRS' [package "sp"] with 1 slot
##   .. .. .. .. .. ..@ projargs: chr NA
##   .. .. .. ..@ srs     : chr "+proj=longlat +datum=WGS84 +no_defs"
##   .. .. .. ..@ history : list()
##   .. .. .. ..@ z       : list()
##   .. ..$ :Formal class 'RasterLayer' [package "raster"] with 13 slots
##   .. .. .. ..@ file    :Formal class '.RasterFile' [package "raster"] with 13 slots
##   .. .. .. .. .. ..@ name        : chr "/Users/jblois/Documents/GitHub/biodata_shortcourse/development/paleoclimate/chelsa_LGM_v1_2B_r2_5m/bio_14.tif"
##   .. .. .. .. .. ..@ datanotation: chr "INT2U"
##   .. .. .. .. .. ..@ byteorder   : chr "little"
##   .. .. .. .. .. ..@ nodatavalue : num -Inf
##   .. .. .. .. .. ..@ NAchanged   : logi FALSE
##   .. .. .. .. .. ..@ nbands      : int 1
##   .. .. .. .. .. ..@ bandorder   : chr "BIL"
##   .. .. .. .. .. ..@ offset      : int 0
##   .. .. .. .. .. ..@ toptobottom : logi TRUE
##   .. .. .. .. .. ..@ blockrows   : Named int 1
##   .. .. .. .. .. .. ..- attr(*, "names")= chr "rows"
##   .. .. .. .. .. ..@ blockcols   : Named int 8640
##   .. .. .. .. .. .. ..- attr(*, "names")= chr "cols"
##   .. .. .. .. .. ..@ driver      : chr "gdal"
##   .. .. .. .. .. ..@ open        : logi FALSE
##   .. .. .. ..@ data    :Formal class '.SingleLayerData' [package "raster"] with 13 slots
##   .. .. .. .. .. ..@ values    : logi(0) 
##   .. .. .. .. .. ..@ offset    : Named num 0
##   .. .. .. .. .. .. ..- attr(*, "names")= chr "offset"
##   .. .. .. .. .. ..@ gain      : Named num 1
##   .. .. .. .. .. .. ..- attr(*, "names")= chr "scale"
##   .. .. .. .. .. ..@ inmemory  : logi FALSE
##   .. .. .. .. .. ..@ fromdisk  : logi TRUE
##   .. .. .. .. .. ..@ isfactor  : logi FALSE
##   .. .. .. .. .. ..@ attributes: list()
##   .. .. .. .. .. ..@ haveminmax: logi TRUE
##   .. .. .. .. .. ..@ min       : num 0
##   .. .. .. .. .. ..@ max       : num 661
##   .. .. .. .. .. ..@ band      : int 1
##   .. .. .. .. .. ..@ unit      : chr ""
##   .. .. .. .. .. ..@ names     : chr "bio14"
##   .. .. .. ..@ legend  :Formal class '.RasterLegend' [package "raster"] with 5 slots
##   .. .. .. .. .. ..@ type      : chr(0) 
##   .. .. .. .. .. ..@ values    : logi(0) 
##   .. .. .. .. .. ..@ color     : logi(0) 
##   .. .. .. .. .. ..@ names     : logi(0) 
##   .. .. .. .. .. ..@ colortable: logi(0) 
##   .. .. .. ..@ title   : chr(0) 
##   .. .. .. ..@ extent  :Formal class 'Extent' [package "raster"] with 4 slots
##   .. .. .. .. .. ..@ xmin: num -180
##   .. .. .. .. .. ..@ xmax: num 180
##   .. .. .. .. .. ..@ ymin: num -90
##   .. .. .. .. .. ..@ ymax: num 90
##   .. .. .. ..@ rotated : logi FALSE
##   .. .. .. ..@ rotation:Formal class '.Rotation' [package "raster"] with 2 slots
##   .. .. .. .. .. ..@ geotrans: num(0) 
##   .. .. .. .. .. ..@ transfun:function ()  
##   .. .. .. ..@ ncols   : int 8640
##   .. .. .. ..@ nrows   : int 4320
##   .. .. .. ..@ crs     :Formal class 'CRS' [package "sp"] with 1 slot
##   .. .. .. .. .. ..@ projargs: chr NA
##   .. .. .. ..@ srs     : chr "+proj=longlat +datum=WGS84 +no_defs"
##   .. .. .. ..@ history : list()
##   .. .. .. ..@ z       : list()
##   .. ..$ :Formal class 'RasterLayer' [package "raster"] with 13 slots
##   .. .. .. ..@ file    :Formal class '.RasterFile' [package "raster"] with 13 slots
##   .. .. .. .. .. ..@ name        : chr "/Users/jblois/Documents/GitHub/biodata_shortcourse/development/paleoclimate/chelsa_LGM_v1_2B_r2_5m/bio_15.tif"
##   .. .. .. .. .. ..@ datanotation: chr "INT2U"
##   .. .. .. .. .. ..@ byteorder   : chr "little"
##   .. .. .. .. .. ..@ nodatavalue : num -Inf
##   .. .. .. .. .. ..@ NAchanged   : logi FALSE
##   .. .. .. .. .. ..@ nbands      : int 1
##   .. .. .. .. .. ..@ bandorder   : chr "BIL"
##   .. .. .. .. .. ..@ offset      : int 0
##   .. .. .. .. .. ..@ toptobottom : logi TRUE
##   .. .. .. .. .. ..@ blockrows   : Named int 1
##   .. .. .. .. .. .. ..- attr(*, "names")= chr "rows"
##   .. .. .. .. .. ..@ blockcols   : Named int 8640
##   .. .. .. .. .. .. ..- attr(*, "names")= chr "cols"
##   .. .. .. .. .. ..@ driver      : chr "gdal"
##   .. .. .. .. .. ..@ open        : logi FALSE
##   .. .. .. ..@ data    :Formal class '.SingleLayerData' [package "raster"] with 13 slots
##   .. .. .. .. .. ..@ values    : logi(0) 
##   .. .. .. .. .. ..@ offset    : Named num 0
##   .. .. .. .. .. .. ..- attr(*, "names")= chr "offset"
##   .. .. .. .. .. ..@ gain      : Named num 1
##   .. .. .. .. .. .. ..- attr(*, "names")= chr "scale"
##   .. .. .. .. .. ..@ inmemory  : logi FALSE
##   .. .. .. .. .. ..@ fromdisk  : logi TRUE
##   .. .. .. .. .. ..@ isfactor  : logi FALSE
##   .. .. .. .. .. ..@ attributes: list()
##   .. .. .. .. .. ..@ haveminmax: logi TRUE
##   .. .. .. .. .. ..@ min       : num 0
##   .. .. .. .. .. ..@ max       : num 331
##   .. .. .. .. .. ..@ band      : int 1
##   .. .. .. .. .. ..@ unit      : chr ""
##   .. .. .. .. .. ..@ names     : chr "bio15"
##   .. .. .. ..@ legend  :Formal class '.RasterLegend' [package "raster"] with 5 slots
##   .. .. .. .. .. ..@ type      : chr(0) 
##   .. .. .. .. .. ..@ values    : logi(0) 
##   .. .. .. .. .. ..@ color     : logi(0) 
##   .. .. .. .. .. ..@ names     : logi(0) 
##   .. .. .. .. .. ..@ colortable: logi(0) 
##   .. .. .. ..@ title   : chr(0) 
##   .. .. .. ..@ extent  :Formal class 'Extent' [package "raster"] with 4 slots
##   .. .. .. .. .. ..@ xmin: num -180
##   .. .. .. .. .. ..@ xmax: num 180
##   .. .. .. .. .. ..@ ymin: num -90
##   .. .. .. .. .. ..@ ymax: num 90
##   .. .. .. ..@ rotated : logi FALSE
##   .. .. .. ..@ rotation:Formal class '.Rotation' [package "raster"] with 2 slots
##   .. .. .. .. .. ..@ geotrans: num(0) 
##   .. .. .. .. .. ..@ transfun:function ()  
##   .. .. .. ..@ ncols   : int 8640
##   .. .. .. ..@ nrows   : int 4320
##   .. .. .. ..@ crs     :Formal class 'CRS' [package "sp"] with 1 slot
##   .. .. .. .. .. ..@ projargs: chr NA
##   .. .. .. ..@ srs     : chr "+proj=longlat +datum=WGS84 +no_defs"
##   .. .. .. ..@ history : list()
##   .. .. .. ..@ z       : list()
##   .. ..$ :Formal class 'RasterLayer' [package "raster"] with 13 slots
##   .. .. .. ..@ file    :Formal class '.RasterFile' [package "raster"] with 13 slots
##   .. .. .. .. .. ..@ name        : chr "/Users/jblois/Documents/GitHub/biodata_shortcourse/development/paleoclimate/chelsa_LGM_v1_2B_r2_5m/bio_16.tif"
##   .. .. .. .. .. ..@ datanotation: chr "INT2U"
##   .. .. .. .. .. ..@ byteorder   : chr "little"
##   .. .. .. .. .. ..@ nodatavalue : num -Inf
##   .. .. .. .. .. ..@ NAchanged   : logi FALSE
##   .. .. .. .. .. ..@ nbands      : int 1
##   .. .. .. .. .. ..@ bandorder   : chr "BIL"
##   .. .. .. .. .. ..@ offset      : int 0
##   .. .. .. .. .. ..@ toptobottom : logi TRUE
##   .. .. .. .. .. ..@ blockrows   : Named int 1
##   .. .. .. .. .. .. ..- attr(*, "names")= chr "rows"
##   .. .. .. .. .. ..@ blockcols   : Named int 8640
##   .. .. .. .. .. .. ..- attr(*, "names")= chr "cols"
##   .. .. .. .. .. ..@ driver      : chr "gdal"
##   .. .. .. .. .. ..@ open        : logi FALSE
##   .. .. .. ..@ data    :Formal class '.SingleLayerData' [package "raster"] with 13 slots
##   .. .. .. .. .. ..@ values    : logi(0) 
##   .. .. .. .. .. ..@ offset    : Named num 0
##   .. .. .. .. .. .. ..- attr(*, "names")= chr "offset"
##   .. .. .. .. .. ..@ gain      : Named num 1
##   .. .. .. .. .. .. ..- attr(*, "names")= chr "scale"
##   .. .. .. .. .. ..@ inmemory  : logi FALSE
##   .. .. .. .. .. ..@ fromdisk  : logi TRUE
##   .. .. .. .. .. ..@ isfactor  : logi FALSE
##   .. .. .. .. .. ..@ attributes: list()
##   .. .. .. .. .. ..@ haveminmax: logi TRUE
##   .. .. .. .. .. ..@ min       : num 0
##   .. .. .. .. .. ..@ max       : num 3276
##   .. .. .. .. .. ..@ band      : int 1
##   .. .. .. .. .. ..@ unit      : chr ""
##   .. .. .. .. .. ..@ names     : chr "bio16"
##   .. .. .. ..@ legend  :Formal class '.RasterLegend' [package "raster"] with 5 slots
##   .. .. .. .. .. ..@ type      : chr(0) 
##   .. .. .. .. .. ..@ values    : logi(0) 
##   .. .. .. .. .. ..@ color     : logi(0) 
##   .. .. .. .. .. ..@ names     : logi(0) 
##   .. .. .. .. .. ..@ colortable: logi(0) 
##   .. .. .. ..@ title   : chr(0) 
##   .. .. .. ..@ extent  :Formal class 'Extent' [package "raster"] with 4 slots
##   .. .. .. .. .. ..@ xmin: num -180
##   .. .. .. .. .. ..@ xmax: num 180
##   .. .. .. .. .. ..@ ymin: num -90
##   .. .. .. .. .. ..@ ymax: num 90
##   .. .. .. ..@ rotated : logi FALSE
##   .. .. .. ..@ rotation:Formal class '.Rotation' [package "raster"] with 2 slots
##   .. .. .. .. .. ..@ geotrans: num(0) 
##   .. .. .. .. .. ..@ transfun:function ()  
##   .. .. .. ..@ ncols   : int 8640
##   .. .. .. ..@ nrows   : int 4320
##   .. .. .. ..@ crs     :Formal class 'CRS' [package "sp"] with 1 slot
##   .. .. .. .. .. ..@ projargs: chr NA
##   .. .. .. ..@ srs     : chr "+proj=longlat +datum=WGS84 +no_defs"
##   .. .. .. ..@ history : list()
##   .. .. .. ..@ z       : list()
##   .. ..$ :Formal class 'RasterLayer' [package "raster"] with 13 slots
##   .. .. .. ..@ file    :Formal class '.RasterFile' [package "raster"] with 13 slots
##   .. .. .. .. .. ..@ name        : chr "/Users/jblois/Documents/GitHub/biodata_shortcourse/development/paleoclimate/chelsa_LGM_v1_2B_r2_5m/bio_17.tif"
##   .. .. .. .. .. ..@ datanotation: chr "INT2U"
##   .. .. .. .. .. ..@ byteorder   : chr "little"
##   .. .. .. .. .. ..@ nodatavalue : num -Inf
##   .. .. .. .. .. ..@ NAchanged   : logi FALSE
##   .. .. .. .. .. ..@ nbands      : int 1
##   .. .. .. .. .. ..@ bandorder   : chr "BIL"
##   .. .. .. .. .. ..@ offset      : int 0
##   .. .. .. .. .. ..@ toptobottom : logi TRUE
##   .. .. .. .. .. ..@ blockrows   : Named int 1
##   .. .. .. .. .. .. ..- attr(*, "names")= chr "rows"
##   .. .. .. .. .. ..@ blockcols   : Named int 8640
##   .. .. .. .. .. .. ..- attr(*, "names")= chr "cols"
##   .. .. .. .. .. ..@ driver      : chr "gdal"
##   .. .. .. .. .. ..@ open        : logi FALSE
##   .. .. .. ..@ data    :Formal class '.SingleLayerData' [package "raster"] with 13 slots
##   .. .. .. .. .. ..@ values    : logi(0) 
##   .. .. .. .. .. ..@ offset    : Named num 0
##   .. .. .. .. .. .. ..- attr(*, "names")= chr "offset"
##   .. .. .. .. .. ..@ gain      : Named num 1
##   .. .. .. .. .. .. ..- attr(*, "names")= chr "scale"
##   .. .. .. .. .. ..@ inmemory  : logi FALSE
##   .. .. .. .. .. ..@ fromdisk  : logi TRUE
##   .. .. .. .. .. ..@ isfactor  : logi FALSE
##   .. .. .. .. .. ..@ attributes: list()
##   .. .. .. .. .. ..@ haveminmax: logi TRUE
##   .. .. .. .. .. ..@ min       : num 0
##   .. .. .. .. .. ..@ max       : num 2029
##   .. .. .. .. .. ..@ band      : int 1
##   .. .. .. .. .. ..@ unit      : chr ""
##   .. .. .. .. .. ..@ names     : chr "bio17"
##   .. .. .. ..@ legend  :Formal class '.RasterLegend' [package "raster"] with 5 slots
##   .. .. .. .. .. ..@ type      : chr(0) 
##   .. .. .. .. .. ..@ values    : logi(0) 
##   .. .. .. .. .. ..@ color     : logi(0) 
##   .. .. .. .. .. ..@ names     : logi(0) 
##   .. .. .. .. .. ..@ colortable: logi(0) 
##   .. .. .. ..@ title   : chr(0) 
##   .. .. .. ..@ extent  :Formal class 'Extent' [package "raster"] with 4 slots
##   .. .. .. .. .. ..@ xmin: num -180
##   .. .. .. .. .. ..@ xmax: num 180
##   .. .. .. .. .. ..@ ymin: num -90
##   .. .. .. .. .. ..@ ymax: num 90
##   .. .. .. ..@ rotated : logi FALSE
##   .. .. .. ..@ rotation:Formal class '.Rotation' [package "raster"] with 2 slots
##   .. .. .. .. .. ..@ geotrans: num(0) 
##   .. .. .. .. .. ..@ transfun:function ()  
##   .. .. .. ..@ ncols   : int 8640
##   .. .. .. ..@ nrows   : int 4320
##   .. .. .. ..@ crs     :Formal class 'CRS' [package "sp"] with 1 slot
##   .. .. .. .. .. ..@ projargs: chr NA
##   .. .. .. ..@ srs     : chr "+proj=longlat +datum=WGS84 +no_defs"
##   .. .. .. ..@ history : list()
##   .. .. .. ..@ z       : list()
##   .. ..$ :Formal class 'RasterLayer' [package "raster"] with 13 slots
##   .. .. .. ..@ file    :Formal class '.RasterFile' [package "raster"] with 13 slots
##   .. .. .. .. .. ..@ name        : chr "/Users/jblois/Documents/GitHub/biodata_shortcourse/development/paleoclimate/chelsa_LGM_v1_2B_r2_5m/bio_18.tif"
##   .. .. .. .. .. ..@ datanotation: chr "INT2U"
##   .. .. .. .. .. ..@ byteorder   : chr "little"
##   .. .. .. .. .. ..@ nodatavalue : num -Inf
##   .. .. .. .. .. ..@ NAchanged   : logi FALSE
##   .. .. .. .. .. ..@ nbands      : int 1
##   .. .. .. .. .. ..@ bandorder   : chr "BIL"
##   .. .. .. .. .. ..@ offset      : int 0
##   .. .. .. .. .. ..@ toptobottom : logi TRUE
##   .. .. .. .. .. ..@ blockrows   : Named int 1
##   .. .. .. .. .. .. ..- attr(*, "names")= chr "rows"
##   .. .. .. .. .. ..@ blockcols   : Named int 8640
##   .. .. .. .. .. .. ..- attr(*, "names")= chr "cols"
##   .. .. .. .. .. ..@ driver      : chr "gdal"
##   .. .. .. .. .. ..@ open        : logi FALSE
##   .. .. .. ..@ data    :Formal class '.SingleLayerData' [package "raster"] with 13 slots
##   .. .. .. .. .. ..@ values    : logi(0) 
##   .. .. .. .. .. ..@ offset    : Named num 0
##   .. .. .. .. .. .. ..- attr(*, "names")= chr "offset"
##   .. .. .. .. .. ..@ gain      : Named num 1
##   .. .. .. .. .. .. ..- attr(*, "names")= chr "scale"
##   .. .. .. .. .. ..@ inmemory  : logi FALSE
##   .. .. .. .. .. ..@ fromdisk  : logi TRUE
##   .. .. .. .. .. ..@ isfactor  : logi FALSE
##   .. .. .. .. .. ..@ attributes: list()
##   .. .. .. .. .. ..@ haveminmax: logi TRUE
##   .. .. .. .. .. ..@ min       : num 0
##   .. .. .. .. .. ..@ max       : num 3276
##   .. .. .. .. .. ..@ band      : int 1
##   .. .. .. .. .. ..@ unit      : chr ""
##   .. .. .. .. .. ..@ names     : chr "bio18"
##   .. .. .. ..@ legend  :Formal class '.RasterLegend' [package "raster"] with 5 slots
##   .. .. .. .. .. ..@ type      : chr(0) 
##   .. .. .. .. .. ..@ values    : logi(0) 
##   .. .. .. .. .. ..@ color     : logi(0) 
##   .. .. .. .. .. ..@ names     : logi(0) 
##   .. .. .. .. .. ..@ colortable: logi(0) 
##   .. .. .. ..@ title   : chr(0) 
##   .. .. .. ..@ extent  :Formal class 'Extent' [package "raster"] with 4 slots
##   .. .. .. .. .. ..@ xmin: num -180
##   .. .. .. .. .. ..@ xmax: num 180
##   .. .. .. .. .. ..@ ymin: num -90
##   .. .. .. .. .. ..@ ymax: num 90
##   .. .. .. ..@ rotated : logi FALSE
##   .. .. .. ..@ rotation:Formal class '.Rotation' [package "raster"] with 2 slots
##   .. .. .. .. .. ..@ geotrans: num(0) 
##   .. .. .. .. .. ..@ transfun:function ()  
##   .. .. .. ..@ ncols   : int 8640
##   .. .. .. ..@ nrows   : int 4320
##   .. .. .. ..@ crs     :Formal class 'CRS' [package "sp"] with 1 slot
##   .. .. .. .. .. ..@ projargs: chr NA
##   .. .. .. ..@ srs     : chr "+proj=longlat +datum=WGS84 +no_defs"
##   .. .. .. ..@ history : list()
##   .. .. .. ..@ z       : list()
##   .. ..$ :Formal class 'RasterLayer' [package "raster"] with 13 slots
##   .. .. .. ..@ file    :Formal class '.RasterFile' [package "raster"] with 13 slots
##   .. .. .. .. .. ..@ name        : chr "/Users/jblois/Documents/GitHub/biodata_shortcourse/development/paleoclimate/chelsa_LGM_v1_2B_r2_5m/bio_19.tif"
##   .. .. .. .. .. ..@ datanotation: chr "INT2U"
##   .. .. .. .. .. ..@ byteorder   : chr "little"
##   .. .. .. .. .. ..@ nodatavalue : num -Inf
##   .. .. .. .. .. ..@ NAchanged   : logi FALSE
##   .. .. .. .. .. ..@ nbands      : int 1
##   .. .. .. .. .. ..@ bandorder   : chr "BIL"
##   .. .. .. .. .. ..@ offset      : int 0
##   .. .. .. .. .. ..@ toptobottom : logi TRUE
##   .. .. .. .. .. ..@ blockrows   : Named int 1
##   .. .. .. .. .. .. ..- attr(*, "names")= chr "rows"
##   .. .. .. .. .. ..@ blockcols   : Named int 8640
##   .. .. .. .. .. .. ..- attr(*, "names")= chr "cols"
##   .. .. .. .. .. ..@ driver      : chr "gdal"
##   .. .. .. .. .. ..@ open        : logi FALSE
##   .. .. .. ..@ data    :Formal class '.SingleLayerData' [package "raster"] with 13 slots
##   .. .. .. .. .. ..@ values    : logi(0) 
##   .. .. .. .. .. ..@ offset    : Named num 0
##   .. .. .. .. .. .. ..- attr(*, "names")= chr "offset"
##   .. .. .. .. .. ..@ gain      : Named num 1
##   .. .. .. .. .. .. ..- attr(*, "names")= chr "scale"
##   .. .. .. .. .. ..@ inmemory  : logi FALSE
##   .. .. .. .. .. ..@ fromdisk  : logi TRUE
##   .. .. .. .. .. ..@ isfactor  : logi FALSE
##   .. .. .. .. .. ..@ attributes: list()
##   .. .. .. .. .. ..@ haveminmax: logi TRUE
##   .. .. .. .. .. ..@ min       : num 0
##   .. .. .. .. .. ..@ max       : num 2944
##   .. .. .. .. .. ..@ band      : int 1
##   .. .. .. .. .. ..@ unit      : chr ""
##   .. .. .. .. .. ..@ names     : chr "bio19"
##   .. .. .. ..@ legend  :Formal class '.RasterLegend' [package "raster"] with 5 slots
##   .. .. .. .. .. ..@ type      : chr(0) 
##   .. .. .. .. .. ..@ values    : logi(0) 
##   .. .. .. .. .. ..@ color     : logi(0) 
##   .. .. .. .. .. ..@ names     : logi(0) 
##   .. .. .. .. .. ..@ colortable: logi(0) 
##   .. .. .. ..@ title   : chr(0) 
##   .. .. .. ..@ extent  :Formal class 'Extent' [package "raster"] with 4 slots
##   .. .. .. .. .. ..@ xmin: num -180
##   .. .. .. .. .. ..@ xmax: num 180
##   .. .. .. .. .. ..@ ymin: num -90
##   .. .. .. .. .. ..@ ymax: num 90
##   .. .. .. ..@ rotated : logi FALSE
##   .. .. .. ..@ rotation:Formal class '.Rotation' [package "raster"] with 2 slots
##   .. .. .. .. .. ..@ geotrans: num(0) 
##   .. .. .. .. .. ..@ transfun:function ()  
##   .. .. .. ..@ ncols   : int 8640
##   .. .. .. ..@ nrows   : int 4320
##   .. .. .. ..@ crs     :Formal class 'CRS' [package "sp"] with 1 slot
##   .. .. .. .. .. ..@ projargs: chr NA
##   .. .. .. ..@ srs     : chr "+proj=longlat +datum=WGS84 +no_defs"
##   .. .. .. ..@ history : list()
##   .. .. .. ..@ z       : list()
##   ..@ title   : chr(0) 
##   ..@ extent  :Formal class 'Extent' [package "raster"] with 4 slots
##   .. .. ..@ xmin: num -180
##   .. .. ..@ xmax: num 180
##   .. .. ..@ ymin: num -90
##   .. .. ..@ ymax: num 90
##   ..@ rotated : logi FALSE
##   ..@ rotation:Formal class '.Rotation' [package "raster"] with 2 slots
##   .. .. ..@ geotrans: num(0) 
##   .. .. ..@ transfun:function ()  
##   ..@ ncols   : int 8640
##   ..@ nrows   : int 4320
##   ..@ crs     :Formal class 'CRS' [package "sp"] with 1 slot
##   .. .. ..@ projargs: chr NA
##   ..@ srs     : chr "GEOGCRS[\"unknown\",\n    DATUM[\"World Geodetic System 1984\",\n        ELLIPSOID[\"WGS 84\",6378137,298.25722"| __truncated__
##   ..@ history : list()
##   ..@ z       : list()
plot(lgm)

plot(lgm[[1]])

## Getting longitude and latitude data

Google “UC Merced GPS coordinates” to find that the place where you’re sitting right now is at longitude 120.42 W, latitude 37.36 N. That is:

UCM_lonlat <- data.frame(lon = -120.42, lat = 37.36)
UCM_lonlat #inspect
##       lon   lat
## 1 -120.42 37.36

(Note that since the longitude is West, we have to make it a negative number.)

Extracting the climate variables

Using the extract() function that we used earlier, get the data from each of the four climate layers you just loaded for the location of UC Merced. (Remember that the syntax is extract(rasterstack,location).) Save them to appropriately named variables.

lgm_UCM <- extract(lgm,UCM_lonlat)
midH_UCM <- extract(midH,UCM_lonlat)
modern_UCM <- extract(modern,UCM_lonlat)
future_UCM <- extract(future,UCM_lonlat)

climate_UCM <- data.frame(rbind(lgm_UCM, midH_UCM, modern_UCM, future_UCM))
rownames(climate_UCM) <- c("lgm","midH","modern","future")
climate_UCM
##           bio1   bio2     bio3     bio4  bio5 bio6  bio7     bio8     bio9
## lgm    11.8000 10.700 36.00000 674.4000 29.80  2.0 29.60  3.90000 21.80000
## midH   17.1000 12.700 41.00000 642.7000 35.10  4.3 30.80 11.10000 25.20000
## modern 16.6015 16.157 46.64261 685.4321 36.14  1.5 34.64 10.02933 24.86667
## future 19.4000 16.700 46.80000 738.5000 39.60  3.9 35.70 11.80000 29.20000
##           bio10  bio11 bio12 bio13 bio14   bio15 bio16 bio17 bio18 bio19
## lgm    22.90000  3.900   572   107     1 83.0000   312     4     5   312
## midH   25.50000  9.900   464    90     0 88.0000   248     3     3   245
## modern 24.86667  8.248   348    64     1 85.4223   187     4     4   174
## future 28.20000 10.400   335    67     1 88.5000   186     4     4   178

Plotting temperature by time

Remember how bioclim variable 1 is mean annual temperature? If you want to see how mean annual temperature at our current location has changed in the last 22,000 years, you can do that now. Let’s plot only that first variable from each of the different time periods.

plot(climate_UCM$bio1)

This is a very simple plot of the data. The points aren’t connected, the axes aren’t labeled, there’s no title, and the x-axis is not scaled by time. We can also do a more sophisticated version by assigning specific age values to the points:

date <- c(-22000, -6000, 0, 50)
plot(x = date, y = climate_UCM$bio1,
     type = "o",
     main = "Temperature at UC Merced", 
     xlab = "Years before present", ylab = "Mean annual temp (degrees C)")

What about the other bioclim variables? Let’s look at precipitation too.

plot(x = date, y = climate_UCM$bio12,
     type = "o",
     main = "Precipitation at UC Merced", 
     xlab = "Years before present", ylab = "Annual precipitation (mm)")

This region has gotten continually hotter and drier since the last ice age.


Exercise

Think of another place on Earth that you feel connected to. (The place where you were born? Your favorite vacation spot? Grandma’s house? Somewhere you hope to live someday?) Repeat the same steps for that place that we just did for UC Merced.

  1. Find its longitude and latitude coordinates. (Remember: North or East = positive numbers, West and South = negative numbers.) Assign those coordinates to a new data frame.
  2. Use extract() to get the values of the bioclimate variables at that location for all four points in time.
  3. Plot the mean annual temperature versus time in your new place, and whichever other bioclimate variables you like. Save the plots to your folder and show them to the instructor.

Part 3. Modeling with paleoclimate data

Now that you’ve looked at both paleoclimate data and paleobiological data, we can put them together to predict future species distributions. We’re going to use the linear model because it worked the best for the last species we tested, but for scientific research, we would usually develop several different models to determine which works best for a specific species.

Note: It may be helpful to save/export each of the plots we generate below for use on Day 3 (as jpeg, png, or tiff files). If you do so, make sure to give your files descriptive names so that you know what the plot is showing.

Loading the fossil occurrence data

load("neotoma_lonlat.RData")
load("gbif.RData")
names(neotoma_lonlat)
## [1] "Microtus pennsylvanicus" "Lynx rufus"             
## [3] "Lontra canadensis"       "Marmota flaviventris"   
## [5] "Erethizon dorsata"       "Bison bison"            
## [7] "Blarina brevicauda"      "Lepus californicus"
names(gbif)
## [1] "Microtus pennsylvanicus" "Lynx rufus"             
## [3] "Lontra canadensis"       "Marmota flaviventris"   
## [5] "Erethizon dorsatum"      "Bison bison"            
## [7] "Blarina brevicauda"      "Lepus californicus"

Open up the file Species List for Day 2.xlsx and take a look at the example species. Choose an organism to investigate in this part (we will try to spread out the species so we don’t all chose the same species).

You can choose from one of the following species: 1. Microtus pennsylvanicus (Eastern meadow vole) 2. Lynx rufus (Bobcat) 3. Lontra canadensis (North American river otter) 4. Marmota flaviventris (Yellow-bellied marmot) 5. Erethizon dorsatum (North American porcupine) 6. Bison bison (American bison) 7. Blarina brevicauda (Northern short-tailed shrew) 8. Lepus californicus (black-tailed jackrabbit)

For the example organism, I’m going to use Microtus pennsylvanicus which is number 1 in this list, but you should use a different one. Assign the index (number) of your organism in the list to index.

index <- 1 #replace 1 with your number here

Now make new variables containing the localities for each time interval and unload the big data files you loaded:

species_name <- names(gbif)[index]
data_lgm <- neotoma_lonlat[[index]]$LGM
data_midH <- neotoma_lonlat[[index]]$MidH
data_now <- gbif[[index]][,c("lon","lat")]
rm(gbif,neotoma_lonlat)

Now you have three variables containing the occurrences of your species for the last glacial maximum (~22,000 years BP), mid-Holocene (~6,000 years BP), and modern (within the last 200 years).

You can map them all on top of each other. (Don’t forget, we have to transform the longitude so that it can be plotted on the map.)

xrange <- transform_lon(range(c(data_lgm$lon,data_midH$lon,data_now$lon),na.rm = T)) + c(-5,5)
yrange <- range(c(data_lgm$lat,data_midH$lat,data_now$lat),na.rm = T) + c(-5,5)

require(maps) #load the mapping library
map("world2",col = "grey",
     xlim = xrange,ylim = yrange)
map.axes()
points(transform_lon(data_now$lon),data_now$lat,col = "maroon1",pch = 20,cex = 0.8)
points(transform_lon(data_lgm$lon),data_lgm$lat,col = "blue",pch = 20,cex = 0.8)
points(transform_lon(data_midH$lon),data_midH$lat,col = "darkorchid",pch = 20,cex = 0.8)
#move the legend to a different corner by changing the value of the argument "x"
legend(x = "bottomleft",legend = c("22,000 ybp (LGM)","6,000 ybp (Mid-Holocene)","Present"),pch = 20,col = c("blue","darkorchid","maroon1"), cex=0.75)

In the Microtus pennsylvanicus example, the vole’s range generally shifts northward through each step.

Note: Do you want to save this plot? If so, give it a descriptive filename and save it for future viewing! For example, the plot below could be saved as: [Yourname]_[SpeciesName]_occurences_through_time.png


Question: Judging the map for your species by eye the same way, do you see a shift in its distribution? What seems to be changing?


Predicting species distributions past and future

We can make a linear model of the species distribution just like before. This time we’ll use all three of the time intervals for which we have data (LGM, mid-Holocene, and modern) to make our predictions.

presvals22k <- extract(lgm, data_lgm[,c("long","lat")])
presvals6k <- extract(midH, data_midH[,c("long","lat")])
presvals0 <- extract(modern, data_now[,c("lon","lat")])
presvals <- data.frame(rbind(presvals22k,presvals6k,presvals0)) #put all presence values together into a single data table
dim(presvals)
## [1] 2309   19

Now we generate the background absence points and extract the climate variable values for them:

set.seed(0) #initialize random number generator
backgr  <-  randomPoints(mask = modern, n = 5000, p = presvals) #this step can take a long time
absvals <- rbind(extract(lgm,backgr),
                 extract(midH,backgr),
                 extract(modern,backgr)) #extract climate values at those points
pb <- c(rep(1, nrow(presvals)), rep(0, nrow(absvals))) #make a vector that indicates whether the species is present at that point or not
sdmdata <- data.frame(cbind(pb, rbind(presvals, absvals))) #bind the presence and absence values together into a single data frame
colnames(sdmdata) <- c("pb",names(modern))

Generate the climate suitability model using the “glm” function.

mod <- glm(pb ~ bio1 + bio2 + bio3 + bio4 + bio5 + bio6 + bio7 + bio8 + bio9 + bio10 + bio11 + bio12 + bio13 + bio14 + bio15 + bio16 + bio17 + bio18 + bio19, data = sdmdata)

Again we use the linear model to make estimates of climate suitability at every location on the map, for each of the time periods (LGM, mid-Holocene, today, and the future):

names(midH) <- names(lgm) <- names(future) <- names(modern)

predlgm <- predict(lgm,mod,progress = "text"); predlgm <- rast(predlgm)
##   |                                                                              |                                                                      |   0%  |                                                                              |======                                                                |   8%  |                                                                              |============                                                          |  17%  |                                                                              |==================                                                    |  25%  |                                                                              |=======================                                               |  33%  |                                                                              |=============================                                         |  42%  |                                                                              |===================================                                   |  50%  |                                                                              |=========================================                             |  58%  |                                                                              |===============================================                       |  67%  |                                                                              |====================================================                  |  75%  |                                                                              |==========================================================            |  83%  |                                                                              |================================================================      |  92%  |                                                                              |======================================================================| 100%
## 
predmidH <- predict(midH,mod,progress = "text"); predmidH <- rast(predmidH)
##   |                                                                              |                                                                      |   0%  |                                                                              |======                                                                |   8%  |                                                                              |============                                                          |  17%  |                                                                              |==================                                                    |  25%  |                                                                              |=======================                                               |  33%  |                                                                              |=============================                                         |  42%  |                                                                              |===================================                                   |  50%  |                                                                              |=========================================                             |  58%  |                                                                              |===============================================                       |  67%  |                                                                              |====================================================                  |  75%  |                                                                              |==========================================================            |  83%  |                                                                              |================================================================      |  92%  |                                                                              |======================================================================| 100%
## 
prednow <- predict(modern,mod,progress = "text"); prednow <- rast(prednow)
##   |                                                                              |                                                                      |   0%  |                                                                              |======                                                                |   8%  |                                                                              |============                                                          |  17%  |                                                                              |==================                                                    |  25%  |                                                                              |=======================                                               |  33%  |                                                                              |=============================                                         |  42%  |                                                                              |===================================                                   |  50%  |                                                                              |=========================================                             |  58%  |                                                                              |===============================================                       |  67%  |                                                                              |====================================================                  |  75%  |                                                                              |==========================================================            |  83%  |                                                                              |================================================================      |  92%  |                                                                              |======================================================================| 100%
## 
pred2070 <- predict(future,mod,progress = "text"); pred2070 <- rast(pred2070)
##   |                                                                              |                                                                      |   0%  |                                                                              |======                                                                |   8%  |                                                                              |============                                                          |  17%  |                                                                              |==================                                                    |  25%  |                                                                              |=======================                                               |  33%  |                                                                              |=============================                                         |  42%  |                                                                              |===================================                                   |  50%  |                                                                              |=========================================                             |  58%  |                                                                              |===============================================                       |  67%  |                                                                              |====================================================                  |  75%  |                                                                              |==========================================================            |  83%  |                                                                              |================================================================      |  92%  |                                                                              |======================================================================| 100%
## 

And plot them:

xrange <- range(c(data_lgm$lon,data_midH$lon,data_now$lon),na.rm = T) + c(-5,5)
yrange <- range(c(data_lgm$lat,data_midH$lat,data_now$lat),na.rm = T) + c(-5,5)

par(mfrow = c(2,2),mar = c(1,3,4,1))
plot(predlgm, main = "22,000 ya",xlim = xrange,ylim = yrange)
points(data_lgm$lon,data_lgm$lat,col = "blue",pch = 20,cex = 0.8)
plot(predmidH, main = "6,000 ya",xlim = xrange,ylim = yrange)
points(data_midH$lon,data_midH$lat,col = "darkorchid",pch = 20,cex = 0.8)
plot(prednow, main = "today",xlim = xrange,ylim = yrange)
points(data_now$lon,data_now$lat,col = "maroon1",pch = 20,cex = 0.8)
plot(pred2070, main = "2070 CE",xlim = xrange,ylim = yrange)

Note: do you want to save/export this plot? If so, go ahead and do that right now.


My example species Microtus pennsylvanicus has had its range expand northward, and the climate variable projections suggest that we should expect this trend to continue. Suitable habitats at the LGM are confined to the South, and as the climate warms and the ice sheets disappear, more habitat is suitable in northern Canada. However, in the 2070 prediction, the current southern part of the range is less green, indicating that it will be less suitable habitat in the near future.


Question: What about your species? What has happened to its range since the last ice age ended, and what do you expect to happen to it by 2070? Do you expect it to be one of the winners or the losers of climate change?


Part 4. Saving your data

Save your data so you can load it again tomorrow. This is not straightforward on UC Merced computer lab computers, so please follow ALL of the following steps:

  1. Choose Session>Save Workspace As…
  2. In the popup window, choose Documents from the list on the left side under Quick Access.
  3. Give the file a UNIQUE name with YOUR NAME in it and click Save. For example, a good file name would be: Blois_Day2.RData (but replace “Blois” with your name!)

Your instructor will make sure these files are here for you to load tomorrow morning.